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ABSTRACT 

In  studying  the  earth  with  reflection  seismics,  one  of  the 
major  unknowns  is  the  velocity  structure  of  the  medium.  Tech- 
niques used  to  determine  the  velocity  structure  commonly  involve 
multi-channel  arrays  which  measure  the  spatial  as  well  as  the 
time  structure  of  the  returning  signals.  The  application  of  a 
data  adaptive  technique,  the  Maximum  Likelihood  Method,  to  the 
problem  of  estimating  seismic  velocities  is  described.  The 
peculiar  problems  of  this  application  are  identified  and  inves- 
tigated. The  windowing  of  short  duration  signals  is  shown  to 
be  an  important  consideration,  and  the  statistics  of  the  MLM 
estimator  for  a single  observation  of  the  data  set  are  presented 
The  adaptive  estimator  is  applied  to  an  ideal  covariance  matrix, 
to  simulated  data,  and  to  field  data.  The  results  show  the  MLM 
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velocity/depth  estimator  to  be  a valuable  tool  in  seismic 
analysis,  and  the  windowing  and  statistical  results  should 
have  general  applications  in  a variety  of  fields. 
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Introduction 


This  thesis  considers  the  application  of  data  adaptive 
array  processing  methods  to  the  estimation  of  velocity/depth 
spectra  in  multi-channel  seismic  reflection  data.  The 
adaptive  processing  methods  are  not  new;  the  basic  techniques 
were  developed  more  than  a decade  ago  for  other  applications, 
and  have  been  applied  to  a multitude  of  time  series  and  array 
processing  problems  to  date.  The  intention  of  the  author  in 
undertaking  this  study  is  to  generalize  the  adaptive  methods 
for  application  to  non-plane  wave,  non-homogeneous  array  data, 
and  to  study  the  requirements  and  performance  of  the  estimation 
methods  as  applied  to  velocity/depth  spectra  estimation.  The 
application  appears  successful  and  the  result  is  an  additional 
tool  for  the  geophysicist  in  his  search  for  higher  resolution 
in  studying  the  earth's  structure.  This  thesis  presents  the 
velocity/depth  spectral  estimators  and  compares  the  conventional 
and  adaptive  forms.  The  details  of  their  implementation  are 
considered  and  an  analysis  of  their  statistics  is  presented. 

The  primary  contributions  of  this  work  are  the  implementation 
of  the  adaptive  processor  to  non-stationary  fields,  the 
importance  and  sensitivity  of  the  time  windowing  to  the 
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conflicting  requirements  of  time  and  frequency  resolution, 
and  an  analysis  of  the  statistics  of  the  adaptive  estimator 
for  a singular  covariance  matrix. 

The  concept  of  remotely  determining  seismic  velocities 
has  been  used  for  many  years  (Green,  1938).  Before  the  advent 
of  the  digital  computer,  the  techniques  involved  the  physical 
manipulation  of  plotted  records  and  the  fitting  of  curves  to 
visually  determined  arrival  times.  Along  with  the  digital 
computer  came  the  ability  to  perform  the  velocity  estimates 
using  correlation  techniques  and  the  reality  of  an  entire 
velocity/depth  spectrum.  A sampling  of  this  development 
may  be  found  in  the  literature  in  papers  by  Green  (1938), 
Durbaum  (1954),  Dix  (1955),  LePichon,  et  al  (1968),  and  Taner 
and  Koehler  (1969) . The  velocity/depth  spectrum  as  we  use 
it  may  be  defined  as  an  estimate  of  the  coherent  reflected 
signal  power  received  from  subsurface  reflectors  as  a function 
of  the  depth  (in  travel  time)  of  the  reflector  and  the  seismic 
RMS  velocity  to  the  reflector.  The  amount  of  effort  expended 
in  velocity  analysis  in  seismics  is  justified  by  the  fact  that 
most  of  the  further  processing  or  analysis  of  data  that  is 
commonly  performed  makes  use  of  the  velocity  information.  In 
particular,  common  depth  point  stacking  and  migration  tech- 
niques depend  heavily  on  accurate  velocity  determinations. 
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Unlike  most  types  of  array  processing,  we  are  dealing 
with  a medium  which  has  non-homogeneous  wave  velocities.  In 
order  to  correctly  phase  or  focus  the  array,  we  must  be  able 
to  relate  the  spatial  position  of  the  array  elements  to  phase 
shifts  or  delay  times.  The  velocity/depth  spectrum  provides 
the  information  which  enables  us  to  do  this.  in  addition 
to  its  applications  in  further  processing  of  the  data,  the 
velocity /depth  information  is  used  in  stratigraphic  interpre- 
tation as  an  aid  in  following  layers  and  in  determining  the 
nature  of  the  structure.  An  important  use  in  geophysical 
interpretation  is  in  differentiating  between  overlapping 
primary  returns  from  deeper  layers  and  multiple  reflections 
from  shallow  reflectors  on  continuous  profiling  records.  The 
normally  higher  velocities  of  deeper  strata  make  it  possible 
to  distinguish  the  two  types  of  returns. 

Data  adaptive  processing  methods  have  been  developed 
in  several  areas  which  include  sonar  array  processing  (see 
■abriel,  1976  for  a good  list  of  references),  time  series 
analysis  (Burg  1967,  Lacoss  1971),  speech  processing  (Makhoul 
1975),  and  communication  theory  (Van  Trees  1968,  Makhoul 
1975) . The  development  was  often  simultaneous,  but  approached 
in  different  ways.  It  is  interesting  to  note  that  although 


each  field  has  its  own  literature  and  terminology  for  the 
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methods,  many  of  them  have  been  shown  to  be  equivalent 
(Edelblute,  et  al  1966,  Gabriel  1976,  Pusey  1975) . Generally 
the  adaptive  methods  may  be  classified  as  one  of  two  types, 
which  have  come  to  be  known  as  the  Maximum  Likelihood  Method 
(MLM)  and  the  Maximum  Entropy  Method  (MEM) . The  MLM  is  attri- 
buted to  Capon  (1967) , but  has  been  shown  to  be  equivalent 
to  several  earlier  techniques  applied  to  single  frequencies 
(Edelblute,  et  al.  1966) . The  MEM  includes  autoregression 
analysis,  covariance  extension,  prediction  error  filters, 
innovations  filters,  and  whitening  filters.  The  MEM  techniques 
are  attributed  to  Burg  (1967)  and  Parzen  (1968,  1969).  Pusey 
(1975)  demonstrates  the  equivalence  of  some  of  the  other  MEM 
forms . 

The  method  we  employ  in  our  study  is  the  Maximum  Like- 
lihood Method.  The  MLM  is  applicable  to  non-homogeneous 
fields  with  non-uniform  sampling,  whereas  the  MEM  has  not  yet 
been  generalized  to  cover  these  cases  in  any  reasonable  manner. 

The  application  of  a data  adaptive  estimation  algorithm  to 
velocity/depth  spectra  estimation  was  first  proposed  by 
Baggeroer  (1974)  and  Baggeroer  and  Leverette  (1975) . This 

I 

thesis  is  a continuation  and  extension  of  that  work. 

The  general  concept  behind  data  adaptive  processing 
methods  is  that  the  filter  coefficients  or  window  weighting 

Ij 

IhT.  •- - - , : : A 


I* 


12 

functions  are  determined  from  the  data  on  each  application 
in  order  to  minimize  the  effects  of  noise  fields.  In  order 
to  demonstrate  this  and  to  further  motivate  a study  of  adap- 
tive array  procedures  applied  to  velocity /depth  spectra 
estimation,  we  would  like  to  give  two  examples.  The  first 
is  an  application  of  the  adaptive  algorithm  to  an  array 
receiving  plane  waves.  The  wave  number  spectra  of  a field 
containing  a single  plane  wave  are  given  in  Figures  1.  and  2. 
Figure  1.  is  the  spectrum  as  measured  by  the  conventional 
array  processor,  and  Figure  2.  is  the  spectrum  as  measured 
by  the  data  adaptive  array  processor.  The  remarkable  increase 
in  resolution  is  more  easily  understood  if  we  examine  the 
beam  patterns  of  the  two  array  processors.  The  conventional 
array  beam  pattern  is  given  in  Figure  3.  For  the  same  array 

with  a noise  field  entering  from  various  directions  (k  ) , the 

N 

adaptive  beam  pattern  is  given  in  Figure  4.  By  adapting  to 
the  received  signal  and  noise  field,  the  adaptive  array  is 
able  to  move  its  peak  and  sidelobes  away  from  interfering 
signals.  This  makes  the  adaptive  processor  particularly 
useful  for  sparse  arrays  which  normally  have  large  sidelobe 
structures.  For  the  second  example.  Figures  5.  and  6.  show 
samples  of  velocity/depth  spectra  generated  by  the  two  esti- 


mators from  data  taken  on  Georges  Bank.  Figure  5.  gives  the 


Figure  3. 


Beam  Pattern  of  Conventional  Array. 


Figure  4.  Beam  Patterns  of  Adaptive  Array  With 
Various  Plane  Wave  Inputs. 
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output  from  a conventional  analysis,  and  Figure  6.  gives 
the  output  from  the  data  adaptive  analysis.  In  the  more 
complicated  case  of  estimating  velocity/depth  information 
instead  of  simple  plane  wave  vectors,  the  adaptive  algorithm 
continues  to  exhibit  a higher  resolution  capability. 

We  begin  with  a review  of  the  travel  time  calculations 
and  the  conventional  estimator.  Although  these  may  be  found 
scattered  throughout  the  literature,  their  importance  to  the 
work  that  follows  and  the  relatively  wide  range  of  audience 
we  hope  to  address  justify  a concise  review.  in  Chapter  1 
we  develop  the  necessary  background  for  the  calculation  of 
propagation  travel  times  from  known  information  about  the 
velocity  structure  of  the  earth.  In  Chapter  2 we  describe 
the  inverse  problem  of  determining  seismic  velocities  from 
measured  travel  times,  making  use  of  the  model  developed  in 
Chapter  1.  The  conventional  estimate  is  presented  in  both 
the  time  and  frequency  domains  and  the  Maximum  Likelihood 
Method  velocity/depth  estimator  is  developed.  Chapter  3 
considers  theoretical  resolution  limits  of  the  conventional 
array  in  terms  of  velocity  and  depth.  The  velocity/time 
ambiguity  function  is  considered,  building  from  the  work  done 
by  Kline  (1976) . Chapter  4 considers  the  problems  encountered 
in  applying  the  estimator  to  real  data,  specifically  the 
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windowing  and  averaging  requirements  in  forming  the  covariance 
matrix.  Chapter  5 develops  the  statistics  of  the  different 
forms  of  the  estimators.  Finally,  Chapter  6 presents  the 
experimental  results  and  conclusions.  The  Appendicies  include 
some  of  the  detailed  calculations  used  in  Chapters  4 and  5, 
and  a glossary  of  symbols. 
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Chapter  1 Array  and  Travel  Path  Geometries  and  Travel 
Time  Calculations. 

Introduction 

Before  addressing  the  problem  of  estimation  of  seismic 
velocities,  it  is  helpful  to  review  some  of  the  physical 
properties  of  the  general  seismic  reflection  problem.  In 
this  chapter  we  review  the  general  array  and  signal  path 
geometries  and  develop  the  commonly  used  RMS  velocity  travel 
time  equation.  The  travel  time,  the  time  required  for  a 
signal  to  traverse  a path  from  the  source  to  a reflecting 
interface  and  back  to  a receiver,  is  one  of  the  most  important 
properties  in  the  estimation  of  velocities.  We  calculate 
the  travel  time  as  a function  of  the  source  to  receiver 
distance  for  a particular  depth  of,  and  RMS  velocity  to  a 
reflecting  surf® re.  We  can  then  generate  a pattern  of  delays 
(or,  in  the  frequency  domain,  phase  shifts)  that  allow  us  to 
steer  or  phase  the  array  to  look  for  coherent  returns  as  a 
function  of  velocity  and  depth. 

Travel  time  calculations  can  become  very  complicated 
for  any  but  the  simplest  geological  models,  and  we  find  that 
simplifications  of  the  geological  models  and  approximate 
solutions  are  desirable  and  necessary  for  our  purposes  in 


velocity/depth  spectra  estimation.  The  RMS  travel  time 
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equation  is  a truncated  series  approximation  of  the  travel 
times  to  interfaces  in  a horizontally  homogeneous  layered 
earth  model.  It  is  a particularly  convenient  model  because 
it  has  a closed  form  solution  and  because  it  simplifies  the 
velocity  dependence  of  the  delay  pattern  to  a single  average 
velocity  rather  than  the  entire  velocity  structure  of  the 
travel  path. 

In  the  remote  measurement  of  seismic  velocities,  we 
measure  the  delay  and  curvature  of  a wavefront  that  has 
originated  at  a point  source  at  the  surface  and  has  penetrated 
the  earth  to  reflect  from  some  lateral  inhomogeneity  in  the 
substrate.  The  most  common  instrumentation  used  to  measure 
the  curvature  of  the  wavefront  is  an  array  of  hydrophones 
or  geophones  uniformly  spaced  along  the  surface  at  increasing 
distances  from  the  source.  The  source  generally  gives  an 
impulsive  signal,  although  longer  coded  signals  which  can 
later  be  deconvolved  or  matched  filtered  are  sometimes  used 
(i.e.,  a chirped  signal).  For  a single  homogeneous  layer 
the  geometry  is  shown  in  Figure  1.1a.  This  is  the  exact 
geometry  for  the  first  return  in  the  case  of  a homogeneous 
and  horizontal  first  layer.  In  marine  data,  it  is  the  water 
column  return  when  there  is  a flat  bottom.  If  the  reflected 
image  is  unfolded  (Figure  1.1b)  and  projected  to  an  array 
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below,  it  is  easily  seen  that  the  wavefront  is  spherical 
and  the  raypaths  are  straight  lines.  The  travel  time  to  a 
receiver  may  be  written  as 


where  C is  the  wave  group  velocity.  We  note  that  for  conven- 
ience and  in  order  to  maintain  consistency,  we  will  use  the 
unit  of  vertical  two-way  travel  time  TQ  to  specify  the 
depth  of  a reflector  throughout  the  remainder  of  this  study. 
Since  the  data  is  a function  of  time,  this  parameter  is  much 
easier  to  correlate  with  the  data  than  would  be  depth  in 
linear  dimensions. 

In  the  case  of  a non-constant  sound  velocity  with  depth, 
we  can  no  longer  assume  straight  line  travel  paths  or  perfectly 
spherical  spreading.  The  rays  will  instead  follow  minimum 
travel  time  paths  as  given  by  Fermat's  Principle.  We  can 
use  Snell's  law  and  ray  path  theory  to  solve  for  the  travel 
time  exactly,  but  the  expression  is  a function  of  the  initial 
angle  and  must  be  solved  parametrically. 

In  order  to  generalize  this  exact  form  of  travel  time 
calculation,  we  consider  a layered  earth  structure  consisting 
of  horizontal  homogeneous  layers.  In  the  limit  as  the  number 
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of  layers  goes  to  infinity  and  the  layer  thicknesses  go  to 
zero,  this  model  may  represent  any  horizontally  homogeneous 
velocity  structure.  The  multi-layer  case  is  depicted  in 
Figure  1.2.  The  ray  parameter  A=  c^/cos (p^  is  preserved  as 
the  wave  travels  through  the  layers.  The  time  through  a 
particular  layer  is 


t 


1.2 


where  tg.  is  the  normal  incidence  travel  time  through  the 
i*1*1  layer.  See  Figure  1.3.  Summing  to  the  m^  layer,  we 
obtain  a two-way  travel  time  of 

r = 2 £ a tj/cj  i. 

**=/ 

The  horizontal  distance  traveled  in  passing  through  each 
layer  is 


-X,  = q tj  5 in  <Pi 


1.4 


Summing  this  over  a two-way  trip  through  m layers  gives  us 


the  total  horizontal  distance  traveled. 
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Figure  1.3  Details  of  One  Layer  of  a Multi-Layer  Travel 
Path.  (Dimensions  in  Seismic  Travel  Time) . 
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1.5 


Given  the  source  to  receiver  distance  Xj  » we  can  solve 
Equation  1.5  for  A.  Inserting  A-  into  Equation  1.3,  we  can 
then  solve  for  the  travel  time  T j . For  the  special  case 
where  the  velocity  in  all  the  layers  is  the  same,  c^  = c^ 
for  all  i,  the  equations  simplify  to 


r 

x 


h 


r 


1.6a 


1.6b 


where 


T.  = Z it., 


l-i 


1.7 


Solving  to  eliminate  A , we  obtain 


T 


a 


1.8 


This  is  identical  to  our  result  for  the  single  layer  case. 

A much  simpler  solution  was  proposed  by  Dix  (1955), 
which  was  a special  case  of  a general  solution  presented  by 
Durbaum  (1954) . A brief  summary  of  the  solution  may  be  found 
in  the  appendix  of  Taner  and  Koehler  (1969) . We  again  refer 
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to  Figure  1.2  for  the  travel  path  geometry  for  a separated 
source  and  receiver.  Following  Taner  and  Koehler,  we  write 
the  travel  time  T as  an  infinite  series  in  powers  of  X * the 
source  to  receiver  distance. 

rx  * A + aX  * * A3X  «-  •••  1.9 


Solving  for  the  first  two  coefficients,  we  obtain 


a.  > lhm  « r 

A,  - [Af/llctt,  - 


1.10a 


1.10b 


An  approximation  using  the  first  two  terms  of  the  series 
gives  us  an  equation  that  is  very  similar  to  the  expression 
for  the  travel  time  through  a single  homogeneous  layer.  If 
we  define 


where  C is  a time  weighted  Root-Mean- Square  velocity,  we 
obtain  a travel  time  expression  of  the  form 
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T 


T 


1.12 


This  is  the  most  common  travel  time  expression  presently  in 

use.  T is  the  two-way  normal  incidence  travel  time  to  the 
o 

layer  of  interest.  C is  known  as  the  RMS  or  stacking  velocity. 
We  note  that  it  is  not  a true  velocity,  but  is  the  first  order 
term  describing  the  hyperbolic  curvature  of  the  wavefront. 

For  normal  array  lengths  and  for  normally  encountered  seismic 
velocity  variations,  the  accuracy  of  this  approximation  for 
the  model  is  better  than  2%  (Taner  and  Koehler,  1969)  . 

If  it  becomes  necessary  to  go  to  the  next  term  in  the 
series,  the  model  becomes  much  more  complicated.  The  coef- 
ficient for  the  next  term  is 


Am  - 


(gc-q*  - (fg(  f t,c;) 

it 


1.13 


Although  we  can  find  no  physical  quantity  corresponding 
directly  to  this  term,  it  is  a measure  of  the  variation  in 
layer  velocities.  A£  goes  to  zero  for  c^  = c^  for  all  i. 
We  expect  this  term  to  be  the  first  order  variation  from  a 
hyperbolic  wavefront  shape.  A ^ may  be  shown  to  always  be 
less  than  or  equal  to  zero  by  the  Schwartz  inequality. 
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The  assumptions  incorporated  in  the  RMS  travel  time 
model  are  the  horizontal  homogeneity  of  the  velocity  structure, 
and  (for  A2  to  be  small)  an  absence  of  extreme  variations  in 
the  vertical  velocity  structure.  In  addition,  all  of  the 
calculations  we  have  considered  so  far  require  that  the  array 
length  be  small  enough  that  there  is  always  a vertical  com- 
ponent to  the  velocity  vector;  that  the  travel  path  does 
not  include  wholly  refracted  segments.  To  put  it  another 
way,  we  must  always  be  close  enough  to  normal  incidence  so 
that  the  interaction  with  the  lowest  interface  is  strictly 
reflection.  As  the  travel  path  deviates  from  vertical,  the 
approximation  in  the  model  becomes  poorer  and  poorer. 

The  most  common  deviation  from  the  assumptions  of  the 
model  is  that  there  is  usually  some  slope  to  the  structure, 
both  in  the  geology  and  the  velocity.  Solving  for  the  first 
order  correction  to  the  model  for  uniform  sloping  layers, 
we  find  that  the  model  is  fairly  robust  to  small  slopes. 

From  model  studies  and  least  squares  fitting  of  real  data, 

Taner  and  Koehler  (1969)  show  that  the  returns  from  mildly 
dipping  layers  are  still  very  closely  hyperbolic  in  form. 
Solving  for  the  delay  times  about  a common  central  ground 
point  for  the  dipping  single  layer  case  (see  Figure  1.4), 


we  obtain 
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The  dipping  layer  always  flattens  the  travel  time  curve  and 
increases  the  apparent  velocity.  Taner  and  Koehler  (1969) 
extend  this  to  multi-layered  cases.  With  all  other  parameters 
held  constant,  increased  dips  produce  higher  apparent  velocities. 
But,  although  the  apparent  velocities  vary,  it  is  important 
that  it  is  still  possible  to  closely  fit  the  delay  pattern 
with  a hyperbolic  model. 

Finally,  we  note  that  it  is  a simple  process  to  take 
the  velocity  structure  in  RMS  velocities  and  calculate  in- 
terval velocities.  The  interval  velocity  between  interface 
i and  i+1  is  given  by 


Summary 

In  the  RMS  travel  time  model  we  have  a simple  and  efficient 
means  of  calculating  the  travel  time  delays  for  the  multi- 
channel array.  The  model  assumes  a horizontally  homogeneous 
acoustic  velocity  structure  for  the  travel  paths,  although 
it  appears  to  be  robust  to  small  dips.  It  is  most  accurate 
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near  vertical  incidence  and  for  structures  without  major 
deviations  in  velocity.  The  model  becomes  invalid  as  any 
part  of  the  travel  path  approaches  a refracting  (i.e.  hori- 
zontal) condition.  With  a means  of  relating  velocity  and 
depth  to  parameters  that  are  directly  measurable,  we  can 
now  look  at  the  estimation  procedure. 
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Chapter  2 Estimation  of  the  Velocity  Function. 

Introduction 

In  this  chapter  we  develop  the  concept  of  a velocity/ 
depth  spectrum  and  present  the  mechanics  of  its  estimation. 

The  form  and  general  structure  of  the  data  are  examined  and 
the  estimation  procedure  is  segmented  into  a two  step  oper- 
ation. The  contribution  of  each  step  toward  the  overall 
resolution  is  examined,  and  areas  of  needed  improvement  ident- 
ified. The  first  step,  the  windowing,  is  shown  to  be  a criti- 
cal, although  often  subtle,  part  of  the  estimation  procedure. 

The  second,  a beamforming  or  coherent  power  estimate,  is  the 
operation  to  which  we  intend  to  apply  the  adaptive  procedure. 

The  conventional  velocity/depth  estimator  is  developed  using 
a beamformer  approach,  and  then  an  adaptive  form  of  this  is 
derived  from  an  adaptive  wave  number  estimator.  Finally,  the 
adaptive  form  is  shown  to  be  computationally  similar  to  the 
conventional  estimator,  and  the  possible  advantage  of  applying 
either  form  in  the  frequency  domain  is  indicated. 

The  Velocity/Depth  Spectrum 

The  concept  of  a velocity/depth  spectrum  has  been  well 
described  in  the  literature  by  LePichon,  Ewing,  and  Houtz  (1968), 
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Taner  and  Koehler  (1969),  and  others.  It  is  an  estimate  of 
the  coherent  power  received  from  a reflecting  surface  at  a 
given  depth  and  at  a given  RMS  velocity.  The  data  set,  com- 
posed of  N channels  of  recordings  from  the  N surface  positions, 
is  scanned  in  an  iterative  process  with  the  estimator.  For 
each  combination  of  depth  and  velocity  the  data  is  windowed 
according  to  the  travel  time  model,  and  an  estimate  of  the 
coherent  power  in  the  windows  is  made  to  form  the  spectral 
estimate.  A sample  spectrum  is  given  in  Figure  2.1. 

There  are  several  ways  commonly  used  to  display  velocity/ 
depth  spectra;  this  one  shows  the  estimated  power  as  the 
displacement  of  plotted  traces.  in  most  of  the  work  which 
follows  we  prefer  to  display  the  spectra  in  contour  plots  of 
the  power  levels  in  6 dB  increments.  Because  of  the  simpli- 
city of  the  equations  and  the  ease  of  correlating  the  spectra 
with  the  original  time  traces  of  the  data,  we  always  consider 
depth  in  the  units  of  seconds  of  two-way  travel  time.  Our 
units  of  velocity  are  RMS  meters  per  second. 

An  idealized  example  of  velocity/depth  spectra  estimation 
is  given  in  Figures  2.2  and  2.3.  Figure  2.2  gives  the  time 
traces  from  8 channels  showing  reflected  returns  from  four 
interfaces.  As  the  data  is  scanned  with  the  estimator,  the 
windows  are  delayed  according  to  a travel  time  model  such 
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Fiqure  2.2  Simulated  Data  Set  Showing  Windows  Properly 
Delayed  for  Third  Reflector. 


Channel  Offset  (meters) 
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Velocity/Depth  Spectral  Estimates  of  Data  in  Fig. 


Travel  Time  (seconds) 
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as  we  calculated  in  the  previous  chapter.  The  windows  in 
Figure  2.2  are  shown  delayed  for  a velocity  and  depth  cor- 
responding to  the  third  reflector.  As  the  velocity  in  the 
travel  time  model  is  incremented  in  the  scanning  process,  the 
window  delays  are  shifted  appropriately.  Examples  of  the  re- 
sulting windowed  data  for  several  shifts  in  velocity  are 
given  in  Figure  2.3a  through  2.3c.  Changes  in  the  depth 
(normal  incidence  travel  time)  shift  the  windows  in  a similar 
manner,  although  much  more  uniformly  up  or  down  the  trace  for 
all  the  channels.  For  each  delay  pattern  specified  by  the 
combination  of  each  depth  and  each  velocity,  the  data  is 
windowed  and  an  estimate  of  the  coherent  energy  in  those  win- 
dows is  made.  The  signals  (though  not  necessarily  the  noise) 
in  the  windows  in  2.3b  are  coherent  across  all  8 channels,  and 
our  estimate  of  the  coherent  power  in  these  windows  will  be 
much  larger  than  the  estimate  for  the  windows  in  Figures  2.3a 
and  2.3c.  This  estimate  of  the  coherent  power  as  a function 
of  the  velocity  and  depth  of  the  delay  model  forms  the 
velocity /depth  spectrum.  The  results  of  the  velocity/depth 
estimation  procedure  for  the  idealized  data  in  Figures  2.2 
and  2.3  are  given  in  Figure  2.4.  The  four  reflectors  are  in- 
dicated by  A and  the  estimates  corresponding  to  the  three 
sets  of  windows  in  Figure  2.3  are  indicated  by  ~}“  . 


40 


Seismic  Reflection  Data 

Before  examining  the  estimation  algorithms  in  any  detail, 
we  first  examine  the  form  of  the  data  and  the  source  signa- 
ture. The  entire  estimation  procedure,  and  the  windowing  in 
particular,  are  ultimately  dependent  on  the  expected  form  of 
the  returning  wave  front.  A typical  example  of  data  is  shown 
in  Figure  2.5.  This  is  data  taken  with  WHOI ' s 6 channel 
system  on  Georges  Bank  in  August  1975.  Reflection  wavefronts 
are  indicated  in  the  time  display  by  hyperbolic  patterns  of 
varying  degrees  of  curvature.  Two  of  these  are  indicated  on 
the  figure.  The  velocity  spectrum  of  this  data  was  given  in 
Figures  5 and  6 in  the  Introduction.  The  set  of  returns 
from  an  interface  is  not  always  obvious,  even  to  the  trained 
eye.  They  vary  for  different  interfaces  and,  to  some  extent, 
from  channel  to  channel.  The  characteristics  of  the  reflected 
wave  are  a function  of  the  source  signature  and  the  dispersive 
and  attenuation  characteristics  of  the  travel  path  medium. 

The  characteristics  of  various  seismic  sources  have 
been  studied  and  classified  (Kramer,  et  al.  1968) . The  out- 
going signal  for  our  data  is  a pulse  from  an  array  of  Bolt 
PAR  airguns.  A typical  outgoing  signature  is  given  in  Figure 
2.6.  It  is  a relatively  wideband  signal  of  approximately 


250  to  500  ms.  duration.  The  frequency  power  spectrum  is 
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given  in  Figure  2.7.  The  spectrum  is  quite  peaked  at  the 
natural  compressional  frequencies  of  the  air  descharge  bubble. 
This  signal  undergoes  phase  changes,  dispersion,  and  selective 
attenuation  as  it  travels  through  the  sediment  structure. 

Since  the  travel  paths  for  the  N channels  of  data  differ  in 
length,  and  usually  to  some  extent  in  composition,  there  will 
be  a modification  of  the  signal  as  a function  of  time  (travel 
distance)  that  will  vary  somewhat  from  channel  to  channel. 

To  the  extent  that  the  signal  from  a given  reflector  is  co- 
herent across  the  array,  our  coherent  power  measurement  func- 
tions well.  Any  incoherence  across  the  wavefront  creates 
difficulties  with  its  measurement  which  we  will  address  later 
when  we  are  considering  the  sensitivity  of  the  estimation 
procedure  to  noise  and  signal  incoherence. 

Partitioning  of  the  Estimation  Procedure 

In  this  section  we  look  separately  at  the  two  basic 
operations  making  up  the  estimation  procedure  - the  windowing 
and  the  coherent  power  estimate.  Each  can  be  used  alone  to 
produce  a form  of  spectrum.  Our  reason  for  doing  so  is  two- 
fold. By  examining  each  aspect  separately,  we  can  better  under- 
stand the  whole  and  how  each  part  contributes  to  the  overall 
resolution  and  accuracy  of  the  complete  estimator.  Secondly, 
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our  proposed  processor  differs  from  the  conventional  in  only 
one  of  the  aspects,  the  coherent  power  estimate,  and  an  under- 
standing of  the  role  of  this  part  enables  us  to  place  limits 
on  any  improvements  we  hope  to  achieve.  In  considering  a 
spectral  estimate  without  the  coherent  power  estimate,  we 
replace  that  operation  with  a calculation  of  the  total  power 
that  is  present  in  the  windows.  For  the  case  of  only  using 
the  coherent  power  estimate,  we  lengthen  the  windows  until 
they  include  the  entire  data  trace.  In  this  way  both  forms 
are  still  estimates  of  the  power  in  the  data  as  a function 
of  velocity  and  depth. 

In  Figures  2.8  and  2.9  we  present  the  two  forms  of 
spectra  run  on  idealized  data  containing  four  reflectors. 
Figure  2.8  gives  a contour  plot  of  the  spectrum  which  relies 
solely  on  windowing  for  its  resolution.  The  points  of  inter- 
est are  the  relatively  sharp  delineation  of  the  reflectors 
in  depth,  but  the  rather  poor  delineation  in  velocity.  Figure 
2.9  gives  the  spectrum  of  the  same  reflectors  calculated  using 
only  the  conventional  coherent  power  estimate.  In  this  case 
there  is  poor  resolution  along  a line  which,  as  we  show  in 
the  next  chapter,  is  defined  by 

2 
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Figure  2.8 


Velocity/Depth  Spectrum  calculated  From 
Incoherent  Arrival  Times.  Four  Reflector 
Simulated  Data  With  No  Noise.  Linear 
Contour  Spacing. 
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Figure  2.9  Velocity/Depth  Spectrum  Calculated  From 
Conventional  Phase  Measurement  Without 
Windowing.  Four  Reflector  Simulated 
Data.  Linear  Contour  Spacing. 
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From  the  general  nature  of  the  two  forms  of  spectra  and 
their  order  of  application,  we  observe  that  the  windowing 
provides  most  of  the  resolution  in  the  time  dimension,  and 
the  coherence  measurement  then  provides  the  resolving  power 
in  the  velocity  dimension.  In  both  of  these  forms  of  spectra 
we  note  that  the  resolution  is  significantly  better  at  shal- 
lower depths.  It  is  interesting  that  the  coherence  measure- 
ment alone  completely  determines  the  reflector  parameters  in 
the  shallowest  region.  The  wavefront  exhibits  the  most 
curvature  (as  determined  by  the  travel  time  equation)  in  the 
very  near  field  of  the  array  and  the  wavefront  shape  is  unique 
for  a given  depth.  In  this  region  the  focusing  of  the  array 
is  analogous  to  holographic  methods.  If  the  entire  geologic 
region  of  interest  were  in  this  holographic  focusing  region, 
we  could  dispense  with  some  of  the  stringent  windowing  require- 
ments. But  such  is  not  often  the  case,  and  we  recall  that 
this  is  also  a region  where  the  travel  time  equations  start 
to  break  down  due  to  refraction  effects.  The  area  where 
we  have  the  most  to  gain  from  new  coherence  measurement  tech- 
niques is  in  the  velocity  resolution  in  the  intermediate  and 
far  end  of  the  Fresnel  region.^  In  these  regions  the  change 

^We  define  the  Fresnel  region  as  being  the  region  where  the 
reflectors  are  shallow  enough  that  the  curvature  of  the 
wavefronts  is  still  significant  over  the  array  length,  and 
the  planewave  approximations  of  the  wavefront  are  not  valid. 
The  term  is  commonly  used  in  optics. 
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in  wavefront  curvature  for  a given  change  in  velocity  is 
relatively  small,  and  any  improvements  in  resolving  power 
effectively  improve  the  resolution  and  the  operating  range 
of  the  array. 


Conventional  Estimator 

In  conventional  array  theory,  a processor  which  calcu- 
lates the  coherent  power  received  by  an  array  is  called  a 
beamformer.  A simple  beamformer  corrects  the  phase  of  the 
signal  from  each  element  to  correctly  "steer"  the  array,  and 
then  sums  the  outputs.  Since  the  phasing  is  a function  of 
frequency,  it  is  often  convenient  to  work  in  the  frequency 
domain.  The  conventional  estimate  of  the  total  coherent 
power  is  given  by 
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where  ^ (V  is  the  frequency  domain  representation  of  the 
signal  from  channel  i, 

and  X IT  n is  the  phase  correction  at  frequency  f for 

channel  i. 

This  estimator  can  be  modified  by  multiplying  each  channel 
by  a weighting  coefficient  in  order  to  taper  the  array,  and 
thus  modify  its  resolution  and  sidelobe  structure.  But  for 
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any  form  of  the  conventional  beamformer,  we  note  that  the 
weights,  and  hence  the  resolution  and  beampatterns,  are 
constant  with  respect  to  the  data  being  looked  at. 

In  the  development  of  velocity/depth  estimation,  the 
traditional  approach  has  been  to  use  an  algorithm  in  the 
time  domain.  We  can  easily  show  that  our  simple  beamformer 
is  equivalent  to  an  un-normalized  "semblance  criteria"  as 
developed  by  Taner  and  Koehler  (1969).  Applying  Parseval's 
theorem  to  Equation  2.1,  we  obtain 

pc  - l 1 1 y,(t*v 

C t H 

The  phase  shifts  become  delays  in  time,  and  the  summation 
in  time  is  over  the  data  window  used  by  the  Fourier  transform 
when  going  to  the  frequency  domain. 

Returning  to  our  frequency  domain  representation,  we 
now  introduce  a vector  notation.  We  let  be  a vector  of 

the  data  ~^(f)  and  E(f)  be  a steering  vector  of  phase  shifts 

j airflj 

6 . Using  this  notation,  the  conventional  estimator 

becomes 

Pc  - £ I 
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The  quantity  [YY*]  is  a matrix  of  products  and  cross  pro- 
ducts of  the  frequency  terms  from  the  Fourier  transforms. 
For  Gaussian  data,  this  is  an  estimate  of  the  covariance 
matrix  of  the  process  (Anderson,  1958) . We  denote  the  co- 
variance  matrix  by  R. 


Rtf)  = Y (0  Y(f) 


2.4 


We  note  that  R(f)  is  hermitian;  it  is  conjugate  symmetric 
complex,  and  is  different  for  each  frequency  of  the  trans- 
form. Collectively,  the  set  of  covariance  matrices  contain 
all  the  relative  phase  information  of  the  N data  windows. 

In  final  form,  we  can  write 


Pc  = E £ 
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Adaptive  Estimator 

The  simple  beamformer  has  a beam  pattern  which  is  directed 
to  look  at  the  amount  of  coherent  energy  in  the  desired  in- 
coming wave  through  the  use  of  the  proper  delays.  The  weights 
on  the  elements  in  this  beamforming  process  are  held  constant, 
so  that  the  basic  shape  of  the  beam  pattern  and  the  associated 
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sidelobe  pattern  for  a given  focus  (velocity  and  depth)  do 
not  change.  But  more  importantly,  they  do  not  depend  upon 
the  data  in  any  direct  manner.  In  order  to  optimize  the 
signal-to-noise  ratio  when  there  are  other  wavefronts  in 
the  viewing  field,  we  would  like  the  beam  pattern  to  adapt 
to  the  data  being  processed.  By  changing  the  weights  of  the 
array  elements,  the  beam  pattern  may  be  controlled  such  that 
the  peak  and  sidelobes  of  the  pattern  are  kept  away  from  the 
directions  that  may  interfere  with  the  estimation  at  a 
particular  desired  direction. 

The  data  adaptive  algorithm  we  are  incorporating  is 
called  the  high  resolution  Maximum  Likelihood  Method,  or  MLM. 
It  was  developed  for  wave-vector  analysis  for  the  large  aper- 
ture seismic  array  (LASA)  in  Montana  by  Capon  (1967) . Our 
application  differs  from  previous  uses  in  that  the  field 
being  measured  does  not  consist  of  plane  waves.  The  data 
field  is  non-homogeneous,  or  spatially  non-stat ionary . This 
characteristic  rules  out  most  other  data  adaptive  methods 
that  are  in  popular  use. 

The  MLM  is  based  upon  the  design  of  a minimum  noise 
unbiased  estimator.  The  estimator  is  constrained  to  pass  the 
desired  wave  (phase  or  delay  pattern)  with  no  distortion. 


while  optimally  suppressing  any  noise  fields.  The  resulting 
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estimator  is  identical  to  the  maximum  likelihood  estimate 
if  the  input  signal  field  is  a multi-dimensional  Gaussian 
process.'*'  The  concept  of  the  MLM  of  wavenumber  estimation 
is  to  calculate  the  average  power  that  this  unbiased,  or 
maximum  likelihood,  estimator  has  as  a function  of  the  steer- 
ing wavenumber,  k.  There  are  several  ways  to  arrive  at  the 
MLM  wavenumber  estimator  formula;  and  we  present  one  which 
has  an  intuitive  appeal  based  upon  the  unbiased  array  pro- 
cessor. Similar  discussion  can  be  found  in  Edelblute,  et  al. 
(1967),  Capon  (1969),  and  Lacoss  (1971). 

The  unbiased  estimator  for  a plane  wave  with  a wavenumber 

k operating  in  the  presence  of  a noise  field  with  a spectral 

2 

cross  correlation  matrix,  R,  is  given  by 

AM  , B(0  EM 

E*(b)  B'(0  £(k)  2-6 

where  R„ (f ) is  the  cross  spectra  between  array  elements 
i and  j at  frequency  f,  and 

^The  maximum  likelihood  estimator  is  the  one  which  gives  as 
its  estimate  the  parameter  set  which  has  the  maximum 
probability  of  producing  the  received  signal.  (see  Van 
Trees,  1968) 

o 

We  use  notation  similar  to  Lacoss  (1971) . 
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is  a steering  vector  consisting  of  the  phase  shifts  required 
for  each  array  element.  Now,  if  the  noise  field  is  applied 
to  the  minimum  variance  unbiased  array  processor,  it  passes 
the  component  in  its  steered  direction  without  attenuation 
and  rejects  the  rest  of  the  field  in  the  manner  which  mini- 
mizes the  output  variance.  Ideally,  then,  the  output  vari- 
ance should  indicate  the  intensity  of  the  component  in  the 
steering  direction,  and  this  is  defined  as  the  MLM  wavenumber 
estimator  formula. 

S^(k)  * a-\h)  = A+(v) 

= [ Ef(U  RK)  EflOj  ' 

The  final  step  is  to  employ  an  estimate  of  the  cross  spectral 
correlation  matrix.'*' 

LM  - [£r«§WE(W ]' 

2.8 


^"Capon  and  Goodman  have  derived  formulae  which  specify  the 
fluctuation  introduced  by  using  an  estimate  of  the  cross 
correlation  matrix.  Essentially,  their  results  show  that 
one  loses  N degrees  of  stability  in  the  MLM  formula  when 
one  has  a multi-dimensional  Gaussian  process. 
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The  form  of  the  MLM  estimator  can  be  compared  with  the  more 
conventional  beamformer  estimator, 

LM  = f eVjRWEfcj 

We  observe  that  additional  computation  required  essentially 
consists  of  inverting  the  cross  spectral  matrix,  which  is  a 
minor  computational  load  when  compared  with  that  of  estimating 
the  matrix  and  scanning  across  the  parameter  set. 

In  modifying  the  MLM  adaptive  spectral  estimation  algo- 
rithm for  use  in  estimating  velocity  spectra,  one  major  modi- 
fication is  required,  and  this  is  the  introduction  of  windows.'*' 
For  depths  or  normal  incidence  times  in  excess  of  that  where 
there  is  holographic  resolution  by  the  phasing  across  the 
array,  the  only  way  that  one  can  obtain  resolution  in  depth 
is  to  use  a sequence  of  window  sets  which  are  positioned  as 
a function  of  depth.  Since  the  velocity  also  influences  the 
position  of  the  windows,  especially  at  the  more  distant  ele- 
ments, these  windows  are  positioned  as  a function  of  both  depth 
and  velocity.  The  net  effect  is  that  one  essentially  has  a 
local  estimate  of  the  cross  spectral  matrix  and  a resulting  MLM 

^Almost  all  previous  applications  of  the  MLM  algorithm  have 
implicity  employed  windows;  but  here  their  role  is  more 
important  because  of  the  inhomogeneity  of  the  spatial  process. 
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velocity  spectra  estimate  around  each  window  position. 

The  presence  of  this  windowing  procedure  introduces  a 
tradeoff  which  turns  out  to  be  quite  important  in  estimating 
the  cross  spectral  matrix.  (In  fact,  understanding  the 
presence  of  this  tradeoff  proved  to  be  one  of  the  more  subtle 
issues  of  this  investigation.)  The  conflicting  issues  in 
this  tradeoff  may  be  summarized  as  follows:  Good  depth  reso- 
lution and  suppression  of  interference  from  reflectors  at 
different  depths  requires  multiplication  by  short  duration 
windows  in  the  time  domain.  This,  however,  implies  a smearing 
of  the  data,  especially  the  phase,  across  the  bandwidth  of 
the  window  which  increases  as  the  window  is  shortened.  We 
analyze  this  tradeoff  in  more  detail  in  Chapter  4. 

With  these  comments  on  the  use  of  windows,  we  define 
the  MLM  velocity/depth  spectra  estimate  to  be 


P (ic^f)  = attc-'f)  gd.of) } 


2.10 


where 


R(rc:f)  = i f X(TPf)Yrncq; 


m*/ 


2.11 
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which  is  an  estimate  of  the  covariance  based  upon  transforms, 
Y (7^C:f)  » °f  the  data  within  windows  positioned  around  depth 

A A 

T and  velocity  C;  and  where  E(T0,C:f)  is  a steering  phasing 
vector  in  the  direction  of  the  desired  depth  and  velocity 

A A 

parameters  Tq  and  C.  If  we  compare  the  form  of  the  MLM 
velocity/depth  estimator  to  the  conventional  beamforming 
procedure  based  upon  coherency  measure,  we  observe  that  it 
is  completely  analogous  to  comparing  the  MLM  and  conventional 
wavenumber  estimators. 

Finally,  we  note  that  the  estimator  is  a function  of 
frequency  and  is  applied  to  discrete  frequency  bands  of  the 
Fourier  transform.  The  characteristics  of  seismic  data  are 
such  that  this  partitioning  of  frequency  is  often  desirable. 
Real  reflecting  horizons  are  often  wavelength  selective 
because  of  the  finite  thickness  of  the  impedance  transition 
region.  Maintaining  separate  estimates  over  frequency  not 
only  gives  sharper  resolution  of  this  type  of  reflection,  but 
gives  some  insight  into  the  nature  of  the  reflecting  surface. 


1 


i 
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Chapter  3 Beam  Patterns  and  Ambiguity  Functions. 

Introduction 

In  the  general  introduction  we  presented  the  on-axis 
beam  patterns  of  linear  arrays  looking  at  a single  plane 
wave  in  wavenumber  space.  These  gave  us  some  insight  into 
the  high  resolution  capabilities  of  the  adaptive  array. 

In  this  chapter  we  examine  the  conventional  beam  pattern 
of  an  array  looking  at  hyperbolic  waves  in  velocity-time 
space.  This  will  provide  us  with  a much  better  indication 
of  the  resolution  of  the  beamforming  process  which  we  looked 
at  in  a superficial  manner  in  the  last  chapter. 

The  general  function  we  need  to  define  this  resolution 
is  the  parameter  ambiguity  function.  The  ambiguity  function 
has  been  described  as  the  response  of  a matched  filter  to 
a mis-matched  signal.  In  the  case  of  an  array  processor, 
it  is  the  normalized  response  of  a steered  array  to  waves 
other  than  the  primary  focus.  We  consider  a uniformly  spaced 
linear  array  as  shown  in  Figure  3.1.  The  array  may  be  steered 
to  receive  waves  from  various  directions  by  adding  appropriate 
delays  to  each  element.  For  plane  waves  and  a linear  array, 
the  ambiguity  function  is  given  by 
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<t>(k,e\k.e.)  - i Z e 


j(/Jksin0  - y(K,0.)) 


3.1 


where  the  steering  function 


ty'MJ 


is  given  by 


- i^sin©. 


3.2 


and  m = (N-l)/2,  N = number  of  elements  in  the  array.  By 
modifying  the  form  of  the  exponent  in  Equation  3.1,  it  is 
easily  seen  that  we  can  form  the  ambiguity  function  for  the 
array  response  to  non-plane  waves.  In  the  case  of  wide  angle 
reflections  from  horizontal  layers,  the  wave  may  be  specified 
by  the  RMS  travel  time  model.  We  then  have  the  ambiguity 
function  in  terms  of  velocity  and  depth. 

,c,|7r,cj=  i L e c*  c' 

i*/ 

This  is  the  complex  monochromatic  ambiguity  function.  For 
a case  where  we  had  a signal  that  was  zero  phase,  we  could 
simply  weight  the  complex  monochromatic  ambiguity  function 
by  the  frequency  spectrum  of  the  signal  and  integrate  over 
frequency  to  obtain  a wide  band  ambiguity  function.  Kline 
(1976)  studied  this  wideband  function  and  found  greatly 


increased  resolution  capabilities.  With  a signal  of  unknown 
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phase  characteristics,  the  ambiguity  functions  (and  signals) 
add  incoherently  across  frequency,  and  we  must  resort  to 
integrating  the  absolute  magnitudes  of  ambiguity  and  signal 
over  frequency.  We  continue  to  weight  by  the  frequency 
spectrum  to  account  for  changes  in  signal  strength.  Our 
wideband  ambiguity  function  becomes 


tft.W.)  - kfjf  S(f)  l<P(xc^c,)j 

B = fdfS(f) 


3.4 


The  velocity-time  ambiguity  function  is  not  solvable  in  closed 
form,  and  thus  requires  numerical  solutions  or  approximating 
functions.  Kline  (1976)  derived  approximations  for  the  peak 
shapes  and  peak  widths  of  this  function  for  monochromatic 
and  narrow  band  cases  which  prove  useful  when  optimizing 
parameters  for  beam  width  or  sidelobe  structure. 

Looking  at  the  monochromatic  ambiguity  function,  we 
find  a large  region  of  ambiguity  stretching  along  a line 
defined  by 


T C 
11 


“ T2C2 


2 


constant 


3.5 
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When  focusing  on  (T-^.C-^),  the  array  will  respond  almost 
equally  well  to  any  return  falling  on  the  line  defined  by 
Equation  3.5.  We  note  that  this  ambiguity  is  independent 
of  frequency.  The  half  power  points  as  approximated  by 
Kline  (1976)  are 


1.8190  1 

f L*J 


3.6 


where  L is  the  equivalent  length  of  the  array.  For  a 
eq 

discrete  element  array,  the  equivalent  length  is 


Leq  = (N-l)d  3.7 

Figure  3.2  gives  a contour  plot  of  an  exact  monochromatic 
ambiguity  function  calculated  for  T^  =2.0  seconds,  C-^  = 2000. 
rn/s,  f = 20  Hz,  and  N = 12.  The  wideband  ambiguity  function 
as  applicable  to  our  data  is  the  sum  of  monochromatic  ambi- 
guity functions  at  discrete  frequency  points  obtained  by  the 
fast  Fourier  transform  of  sampled  data.  The  general  form  of 
the  ambiguity  function  is  not  changed,  although  the  peak  is 
better  defined.  An  example  is  given  in  Figure  3.3  for  T^  = 

2.0  seconds,  = 2000.  m/s,  f = 20.,  24.,  28.,  and  32.  Hz, 


and  N=  12 . 


Velocity  (m/s) 


'OOfil 
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Discussion 

We  see  that  the  array  focusing  - the  coherent  power 
estimate  - allows  us  to  resolve  a reflection  return  to  a 
one  dimensional  strip  or  line  in  velocity-time  space.  We 
depend  on  the  time  windowing  to  provide  resolution  along 
the  length  of  this  strip.  The  effect  of  applying  the  adap- 
tive processor  will  be  primarily  to  reduce  the  width  of  the 
strip.  The  time  windowing  will  continue  to  carry  the  load 
of  resolution  along  the  length  of  the  strip.  Going  to  a 
wide  band  estimator  does  not  produce  any  significant  improve- 
ments in  the  ambiguity  function.  Higher  frequencies  give 
improved  resolution,  but  our  primary  reason  for  applying 
a wideband  estimator  will  be  for  improved  signa 1-to-noi se 
ratios.  In  the  next  chapter  we  investigate  the  windowing 
to  remove  the  ambiguity  along  the  strip,  and  in  Chapter  6 
we  see  how  reducing  the  width  of  the  strip  greatly  enhances 
the  overall  resolution  of  the  velocity/depth  estimator. 
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Chapter  4 Estimation  of  the  Cross  Spectral  Correlation 
Matrix . 


Introduction 

Both  the  MLM  and  the  frequency  domain  implementation 
of  the  conventional  semblence  criteria  for  estimating  velo- 
city/depth spectra  involve  determining  the  cross  spectral 
correlation  matrix  in  one  way  or  another.  in  applications 
to  stationary  homogeneous  signal  fields  this  typically  in- 
volves averaging  over  transformed  segments  of  the  data  from 
each  of  the  channels.  in  the  application  to  velocity /depth 
spectra,  however,  the  transient  nature  of  the  reflected 
signals  requires  a windowing  operation,  particularly  for 
resolving  along  the  depth,  or  time  coordinate.  The  details 
of  the  cross  spectral  correlation  matrix  estimation  involving 
this  windowing  operation  are  critical,  for  the  errors  and 
biases  introduced  propagate  directly  into  the  final  spectral 
estimate.  The  estimation  of  this  matrix  has  proven  to  be 
the  most  subtle  aspect  of  our  experiments  in  applying  the 
MLM  to  velocity/depth  spectra  estimation. 

The  procedure  for  estimating  the  cross  spectral  corre- 
lation matrix  using  a window  is  shown  schematically  in  Figure 
4.1.  At  a given  frequency  the  diagonal  components  of  this 
matrix  are  measures  of  the  energy  at  each  channel,  while  the 
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off-diagonal  terms  are  indicative  of  the  coherent  energy 
and  its  relative  phasing  from  channel  to  channel.  The  two 
most  important  aspects  in  the  estimation  of  these  components 
are  the  smearing,  or  bias,  and  the  variance.  As  in  any 
spectral  estimation  problem  there  are  inevitable  tradeoffs 
between  these  two  quantities;  the  windowing,  however,  further 
complicates  this  issue.  In  this  chapter  we  examine  some 
aspects  of  estimating  this  matrix  - both  the  smearing  intro- 
duced by  the  windowing  and  the  various  ways  of  averaging  to 
improve  the  stability  of  it. 

4 -A  Windows  and  the  Bias  of  Transforms 

The  spectral  correlation  matrix  is  estimated  using  the 
direct  or  FFT  method  of  spectral  analysis,  so  the  first  step 
involves  analyzing  the  bias  introduced  by  windowed  Fourier 
transforms.  In  this  section  we  examine  this  by  first  intro- 
ducing a stochastic  model  for  the  reflected  signal  from  which 
we  can  calculate  bias  errors  using  established  methods  of 
spectral  analysis.  Then  we  examine  the  effects  of  windowing 
on  an  airgun  source  signature  which  ideally  should  be  repre- 
sentative of  the  signal  reflected  from  a horizon.  Finally, 
we  use  estimates  of  allowable  positional  errors  determined 
from  the  results  of  the  stochastic  analysis  to  derive  bounds 


I 


68 

on  perturbations  of  the  moveout  parameters  T , C for  main- 
taining a particular  level  of  average  bias  in  the  windowed 
transforms. 

We  model  the  reflected  signal  observed  at  an  array  element 
as  a desired  signal  plus  an  additive  noise,  or 

Y (t)  - A(t-\)  t n(t)  4.1 


where  A,  (t)  is  the  reflected  signal  at  the  array  element 
which  arrives  with  a total  travel  time  delay 
or  moveout  of  7^. 

n(t)  is  an  additive  noise  which  may  include  both 
ambient  noise  and  reverberation  from  other 
horizons. 

As  indicated  in  Figure  4.1,  the  windowed  transform  operation 
consists  of  multiplying  the  signal  by  a window  function 
centered  at  T and  then  Fourier  transforming,  or 


- r 00 

Jjf)  = j Y(t)w(t-Z)e 


-jirrft 


(We  use  continuous  time  notation,  although  in  practice  the 
FFT  algorithm  is  used.)  We  specify  the  windows  to  have  a 
half  width  duration  of  M seconds,  and  some  commonly  employed 
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windows  are  indicated  in  Figure  4.2. 

If  there  were  no  windowing  and  no  noise,  i.e.  w(t)  = 1., 
n(t)  = 0.,  the  result  of  the  transformation  would  be 

a 

(fj  = */(f)  & 4.3 


which  consists  of  the  desired  signal  transform  and  a linear 
phase  shift  from  the  travel  time  delay.  Both  the  windowing 
and  the  additive  noise  term  introduce  errors  in  this,  so  one 
actually  obtains 

- | (Mt-i)  + n(t)j  w(t-Tw)  e dt  4.4 

-00 

It  is  convenient  at  this  point  to  define  the  error,  since 
this  is  what  we  wish  to  quantify.  We  have 

fw  -\7jrft 

E(f)  = I [Mt-Z)(w(t-X ,)-i)  + n(t)w(t-rj]  e dt  4.5 

-CD 

Qualitatively,  the  duration  of  the  window,  M,  introduces  a 
tradeoff.  A long  window  leads  to  low  resolution  of  the 
depth  and  higher  noise  in  the  transformation;  however,  it  is 
relatively  insensitive  to  its  exact  positioning  and  intro- 


duces little  bias  or  smearing.  Conversely,  a short  window 
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Figure  4.2a  Commonly  Used  Windows  in  the  Time  Domain 


Rtctangulor 
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leads  to  higher  resolution  in  depth  and  lower  noise;  however, 
it  is  very  sensitive  to  its  exact  positioning  and  can  intro- 
duce significant  bias,  or  smearing  of  the  frequency  domain 
signal . 

For  our  stochastic  analysis  we  model  the  reflected 
signal  as 

A{t)  = a(t)  At)  4-6 


where  Cl(t)  is  an  envelope  function  of  approximate  duration 
^ (half  width)  which  models  the  transient,  or 
short  duration  nature  of  the  reflected  signal; 
*ft)  is  a wideband  stationary  process  which  models 
the  waveform  variation  of  the  signal  within 
the  duration  of  the  envelope. 

If  we  assure  that  the  signal  and  the  noise  are  uncorrelated 
processes,  we  can  determine  the  mean  square  error  by  squaring 
and  averaging  Eq.  4.5.  If  we  express  all  of  the  correlations 
in  terms  of  their  associated  spectra,  we  obtain 


<»r  = J Sjvj  ja(t-Z)(w(t-z)  -l)  eJ 


-j2ir(f-v)t 

dt 


civ 


4.7 
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w(t -rw)  e it 


dv 
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where  and  on  are  the  power  density  spectra  of  the  processes 
x(t)  and  n(t)  respectively.  We  next  assume  that  these  spectra 
are  essentially  constant  across  the  bandwidths  of  the  window 
w(t)  and  the  envelope  Gf(t).  (This  is  a common  assumption  in 
spectral  analysis.)  We  then  can  take  them  outside  the  integrals, 
and  after  using  Parseval's  theorem  we  obtain 


r 

J-O, 

•+ 


5,(f)  f jo(t-Z)(w(t-V-0)  Jt 
5„(f]J  w*(t)  Jt 


4.8 


The  details  of  this  derivation  are  given  in  Appendix  I.  The 
first  term  describes  the  error  introduced  by  the  duration  and 
position  of  the  window  with  respect  to  the  desired  reflected 
signal,  while  the  second  term  describes  the  effects  intro- 
duced by  the  additive  noise.  We  consider  each  of  them  separ- 
ately. 

The  noise  term  is  easy  to  analyze.  For  almost  any 
reasonable  window,  one  can  demonstrate  that 


oo 

= s„(f)K w/n 

-a? 


4.9 


where  ttw  is  a window  factor  whose  precise  value  depends 
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upon  the  shape  of  the  window,  but  typically  ranges  from 
*5<Kw<2*  The  most  important  observation  is  that  the  RMS 
value  of  the  noise  increases  as  J/v\  , so  one  wants  to  avoid 
excessively  long  windows  for  noise  as  well  as  resolution 
considerations . 

The  signal  term  is  generally  the  more  important  one, 
and  it  is  somewhat  more  difficult  to  analyze.  First  it  is 
convenient  to  normalize  it  simply  for  the  purposes  of  com- 
parison. The  mean  square  value  of  the  desired  signal  with 


no  windowing  is  given  by 


| J(t) 


ss  Sx( f)  ( V( t)  dt 

J-  CO 


4.10 


The  expression  which  quantifies  the  relative  effects  of  the 
mean  square  bias  error  due  to  the  windowing  is  then  given  by 


t\#AX)  - | [a(t)(w(t.r)-lj]W j a\t)Jt 


4.11 


where  at-Xl-X,  is  the  difference  between  the  position  of 
the  window  and  the  center  of  the  desired  signal.  The  precise 
shape  of  this  function  depends  upon  the  particular  window 
and  envelope  employed.  Figures  4.3  and  4.4  are  indicative 
of  the  general  structure.  Figure  4.3  was  computed  using 


Versus  Position  Error  for  Various  Ratios  of 
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i 


Gaussian  shaped  functions  for  the  window  and  envelope  of  the 

4.12a 

4.12b 

Figure  4.4  was  computed  using  Hanning  windows  for  the  shape 
of  both  the  window  and  envelope  of  the  form 

- 1 

w(t)  = i(l+  cos(a?)) 

a{t)  = i (l  + cos(^)) 

Essentially  these  figures  suggest  that  for  less  than  a con- 
servative 10%  error  in  the  average  bias  of  the  windowed 
transform  operation,  one  wants  to  keep  the  positional  error 
within  ±0.1  (i.e.  20%  of  the  effective  window  extent)  and 
use  windows  with  1^/M  ^0.5,  i.e.  windows  whose  duration  is 
at  least  twice  the  effective  signal  duration. 

To  test  the  effects  of  windowing  on  actual  data,  several 
tests  were  performed  upon  recorded  airgun  signatures,  which 


4.13a 


4.13b 


Figure  44  RMS  Error  Versus  Position  Error  for  Various  Ratios 
of  Envelope  Duration  to  Window  Duration 
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ideally  should  represent  the  signal  reflected  from  > horizon. 
The  signature  and  its  unwindowed  transform  are  illustrated 
in  Figure  4.5.  One  can  estimate  the  energy  distribution 
about  a central  location  by  calculating  the  median  signal 
location  and  then  computing  the  residual  energy  outside  an 
interval  about  that  point.  This  suggests  that 

Figure  4.6  illustrates  the  windowed  transform  with  no 
positional  error  using  the  windows  indicated  in  Figure  4.2 
with  a value  of  M = 0.128  secs.,  i.e.  not  conforming  to  our 
previously  suggested  design  guideline  of  t^/M  0.5.  One 
can  observe  that  there  is  some  evident  spectral  smearing, 
but  the  windowed  trasform  is  basically  accurate.  (One  needs 
to  compensate  visually  for  the  phase  jumps  at  ±1Tas  a shift 
of  2ft  in  phase  is  equivalent.)  Figure  4.7  is  a more  sensitive 
indication  of  the  accuracy  of  the  windowed  transform  with 
respect  to  positional  error.  Here  we  have  plotted  the  phase 
deviation  from  linearity  for  the  10  Hz  component  as  the  sig- 
nature is  delayed  through  the  window.  We  can  observe  that 
for  only  two  of  the  windows  is  there  a comparatively  narrow 
range  of  ±0.012  sec.,  otAT/M&O.I  where  the  phase  deviation 
is  within  ±15°  for  an  error  of  *»*30%.  This  is  essentially  in 
agreement  with  Figure  4.3  which  predicts  that  for  'l^/M  — 1., 
the  error  should  be  constant  at  28%  for  AT/t\<0.1,  and  then 
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Figure  4.5a  300  cu.  in.  Airgun  Signature  (including 

water  surface  image) . 


A)  Time  Signature 


Time  (seconds) 


Figure  4.5b  and  4.5c  Frequency  Signature  of  Airgun. 

B)  Frequency  Magnitude 


Frequency  (Hz) 
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Figure  4.6  Frequency  Spectra  Estimates  of  Gun 

Signature  Using  Various  Lag  Windows. 


Frequency  (Hz) 
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Figure  4.7  Variation  From  Linear  Phase  Shift  of 
Gun  Signature  Moving  Through  Various 
Windows . 


Delay  of  Signature  in  Window 


( seconds) 


4 


increase  significantly  thereafter.  Figure  4.4  gives  similar 
results  for  the  Hanning  windows. 

The  final  step  in  our  analysis  of  bias  error  introduced 
by  windowing  ie  to  translate  the  tolerance  in  positional 
error  to  allowable  perturbation  in  normal  incidence  time,  Tq, 
and  velocity,  C.  Essentially,  we  have  that  if  changes  in 
these  parameters  produce  large  positional  arrors  around  a 


normal  moveout  curve  for  the  array  elements,  then  we  require 
a dense  scanning  in  estimating  the  spectral  correlation 
matrix.  Obviously  this  is  an  added  computational  burden 
which  one  would  like  to  avoid. 


We  can  perform  this  analysis  by  taking  the  total  deriv- 
ative of  the  normal  moveout  relationship. 


Ti(j0C}  xj  - /"C*  + (x,-/c): 


4.14a 


ATL  = Tat;  - (xf/G>G 

Ti 


4.14b 


This  can  be  manipulated  into  the  form 


4T(t,C,Xi)  = cos<t>.  *T  - T-«<fcsm<k  X AC  4>15 


where 


T.„'7  x,/cx] 


4.16 
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The  easiest  way  to  employ  this  relation  is  to  note  that  the 
maximum  effects  of  a change  in  AX  are  when  (p.^.0.  and  in 
UC  when  (f){& 90°.  We  can  use  a worst  case  analysis  for  a 
nominal  Tq,  C by  considering  the  situations  at  0- =. 0 . and 
<t>,  - 1,  - iW5T'l  where  Xmax  is  the  array  element 
with  the  most  distant  offset. 

As  a simple  example  we  consider  a velocity  analysis  for 
a 2.5  km.  array  at  Tq  = 3.  secs,  and  "C  ranging  from  1.5 
km/sec  to  4.5  km/sec.  From  our  previous  analysis  we  allow 
a positional  error  of  ±0.025  secs  which  is  divided  equally 
between  that  caused  by  and  that  by  We  then  have 

for  the  allowable  normal  incidence  time  change 

G os  <t>  ! rvtax  All  ^ .012.5"  acts 

or  I wjax  *T0  I ^ .0/2.5"  sac.s  . 4.17 
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, | &E  [ - .0X3  U/se c;  4.19 

while  at  “P  = 4.5  km/sec  it  implies 

^max  * 10’^°  i I max  I ~ .S5i>  kh/sec.  4.20 

Obviously  the  positional  errors  are  more  sensitive  at  the 

lower  velocity,  which  requires  a denser  selection  of  nominal 

parameters  for  T , C.  The  results  for  velocity  increments 
o 

for  the  same  2 . 5 km  array  for  a range  of  velocities  and 
depths  are  plotted  in  Figure  4.8.  Note  the  large  increments 
that  are  allowed  in  the  deep,  high  velocity  region. 

Discussion 

We  now  can  set  up  the  iteration  over  velocity  and  depth 
in  an  optimum  manner.  We  scan  the  estimator  on  increments 
corresponding  to  the  finest  resolution  that  we  can  expect  in 
the  given  dimension.  In  time  the  increment  is  determined 
by  the  length  of  the  signature  and  the  length  of  the  data 
window.  in  velocity  it  is  dependent  on  the  array  spacing 
and  length,  on  the  frequency,  and  on  the  estimator  form  used. 
The  results  of  i-is  last  section  offer  relief,  however,  from 
the  necessity  of  having  to  form  a new  covariance  matrix  for 
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Figure  4.8  Maximum  Velocity  Estimate  Increments  versus 
Velocity  and  Depth  of  the  Estimate  for  a 
2 . 5 km  Array. 
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each  increment  of  the  estimator.  The  estimate  can  be  per- 
formed on  the  same  matrix  without  appreciable  degradation 
over  a range  specified  by  Eqs.  4.16  and  4.17,  and  Figure  4.8. 

4-B  Averaging  and  the  Stability  of  the  Cross-Spectral 

Correlation  Matrix. 

In  the  previous  section  we  concentrated  upon  producing 
an  estimate  of  the  cross  spectral  correlation  matrix  which 
had  a minimum  of  bias.  In  this  section  we  consider  the  other 
aspect  of  this  estimate,  that  of  its  variance  or  stability. 

The  stability  of  the  estimate  is  essentially  determined  by 
the  deterministic  components  and  the  available  number  of  in- 
dependent degrees  of  freedom  in  reducing  any  random  components . 
The  deterministric , or  mean,  components  are  indicative  of 
the  presence  of  reflection  horizons,  while  the  random  ones 
represent  the  variation  that  one  observes  in  the  reflections 
from  them.  The  random  components  may  be  caused  by  variations 
between  the  travel  paths  of  adjacent  shots,  dispersion  be- 
tween different  frequencies,  or  errors  caused  by  random 
noise.  In  this  section  we  examine  the  methods  by  which  one 
may  increase  the  stability  of  the  estimate.  We  reserve  until 
the  following  chapter  a discussion  of  the  statistics  and 
probability  models  for  the  estimators.  The  probabalistic 
models  for  describing  a non-linear  estimator  with 


many  controlling  parameters  tend  to  become  intractable. 

By  first  considering  how  one  goes  about  stabilizing  the 
estimate,  we  gain  some  insight  into  the  description  of  the 
statistics  of  the  complete  estimator  which  we  examine  in 
Chapter  5. 

The  primary  mechanism  for  increasing  the  stability  of 
a spectral  estimate  is  one  of  averaging  over  blocks  of  data. 
Within  the  constraints  of  our  windowing  requirements  there 
are  two  domains  over  which  one  can  average  to  reduce  the 
variance  of  the  estimate  - across  shots  and  across  frequency. 
This  averaging  of  the  data  may  be  performed  at  several  positions 
before,  within,  and  after  the  application  of  the  estimation 
procedure,  each  with  slightly  differing  results.  These 
positions  are  indicated  in  Figure  4.9.  The  two  averaging 
domains  are  sufficiently  different  from  each  other  that  each 
bears  a separate  set  of  comments. 

The  possibility  of  averaging  over  successive  shots  is 
suggested  by  the  similarity  of  signals  produced  by  closely 
spaced  shots.  The  estimate  is  improved  only  if  the  signals 
being  processed  are  coherent  in  some  respect  across  the  shots 
being  averaged,  and  the  noise  is  uncorrelated.  The  effective- 
ness of  this  is  then  a function  of  the  horizontal  homogeneity 

I 

of  the  medium  and  the  distance  between  shot  points,  as  well 
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as  the  control  of  the  array  geometry  and  the  stability  of  the 
airgun  signatures.  For  a horizontal  planar  structure  and 
closely  spaced  shot  points,  the  signal  may  be  coherent  over 
many  shots  and  extensive  averaging  is  possible. 

There  are  two  respects  by  means  of  which  a signal  may 
be  coherent  over  a shot  sequence.  in  the  first  the  wave- 
form may  repeat  from  shot  to  shot,  and  here  a linear  aver- 
aging of  signals,  or  their  transforms  (since  the  transfor- 
mation is  a linear  operation)  is  appropriate.  This  is  indi- 
cated in  the  first  two  averaging  columns  of  Figure  4.9. 
Alternatively,  the  signal  may  very  from  shot  to  shot,  but 
the  correlation  and  relative  phasing  may  be  stable.  Here 
a quadratic  averaging  of  the  cross  products  used  in  estimating 
the  cross  spectral  correlation  matrix  is  appropriate.  This 
is  illustrated  in  the  third  averaging  column  of  the  figure. 

In  estimating  the  matrix  one  can  average  across  fre- 
quency if  the  signals  are  broadband,  and  the  relative  phasing  is 
not  severely  distorted  across  the  frequency  band  used.  The 
same  concepts  that  appear  in  the  analysis  of  conventional 
planar  arrays  also  appear  here.  (See  Skolnik,  1962).  The 
basic  calculation  that  is  performed  is  to  compute  the  bandwidth 
of  the  array  about  a nominal  center  frequency.  For  the  case 
of  a simple  linear  array  in  a field  consisting  of  a single 
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plane  wave,  we  have  a normalized  response  given  by 

y ^i,  t&,)  y Z!  y 0 €.  j 


4.21 


where  the  steering  function  is  given  by 


= Id  ^ 0, 


4.22 


This  is  the  plane  wave  ambiguity  function,  (see  Eq.s  3.1 
and  3.2.)  For  a correctly  steered  array,  we  have 

M. 

C,*C.  4.23 

e,  - e. 

Now,  if  we  let  fQ  vary  while  keeping  f^  fixed,  we  have 

’ctxir  .. 


f(f.,c.e  | C,C„©.)  = jIt  Z e 


4.24 


The  response  to  waves  of  other  frequencies  is  dependent  on 
the  propagation  velocity  and  the  angle  of  incidence,  as  well 
as  on  the  frequency  shift.  For  the  case  of  the  bandwidth  of 
the  array  in  velocity/depth  estimation,  the  response  to  other 
frequencies  is  given  by 


4.25 


M j«rT(-n,e,x,)(f.-() 

f(£.x.C|f,i:.c;  = ip 

Because  of  the  complexity  of  the  geometry  here,  it  is  diffi- 
cult to  state  anything  very  general.  The  array  bandwidth 
depends  on  the  depth  and  velocity  of  focus,  as  well  as  on  the 
array  geometry.  For  the  adaptive  processor,  the  increased 
resolution  will  decrease  the  array  bandwidth  significantly, 
although  this  is  even  more  difficult  to  quantify. 

These  comments,  however,  do  not  hold  for  averaging  a- 
cross  frequency  in  column  5 of  Figure  4.9;  averaging  the  final 
estimate.  Averaging  at  this  point  produces  a wideband  esti- 
mate as  described  in  Chapters  2 and  3.  Here  there  is  no 
longer  phase  information  and  the  estimates  average  coherently 
as  long  as  the  information  in  the  two  frequency  bands  is 
consistent . 

We  have  found  that  in  regions  with  a reasonable  amount 
of  horizontal  homogeneity,  the  velocity/depth  spectra  are 
quite  consistent  across  adjacent  or  closely  spaced  shots. 

(See  Chapter  6)  Averaging  across  shots  in  any  of  the  posi- 
tions is  of  some  benefit.  We  have  had  mixed  results,  however, 
in  averaging  across  the  frequency  domain.  In  all  of  the 
positions  except  column  5 the  smearing  has  been  noticeable. 

In  column  5 we  have  found  that  it  is  often  useful  to  main- 
tain the  separate  estimates  over  frequency.  It  appears  that 
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the  reflection  process  can  be  frequency  selective,  with 
horizons  which  are  evident  in  a CDP  profile  appearing  only 
in  some  of  the  velocity/depth  estimates  versus  frequency. 

It  may  be  possible  to  use  this  frequency  selectivity  con- 
structively, either  for  the  design  of  filters  in  subsequent 
stacking  operations  or  as  a diagnostic  tool  in  interpreting 
the  character  of  the  reflection  horizons. 

Summary 

The  estimation  of  the  covariance  matrix  from  the  data 
is  a critical  step  in  forming  the  velocity/depth  spectrum. 

Two  important  aspects  of  this  estimation  are  the  time  windowing 
prior  to  transforming  and  the  averaging  of  the  data.  The 
optimum  window  shape  and  length  is  dependent  on  the  reflected 
signature.  Once  the  window  is  determined  from  a tradeoff  of 
time  resolution  and  frequency  smearing,  the  bias  due  to 
positional  errors  is  easily  calculated.  Defining  limits 
for  this  bias,  we  can  then  perform  the  velocity /depth  spectral 
power  estimate  over  a small  range  of  depths  and  velocities 
using  the  same  estimate  of  the  covariance  matrix.  Numerically, 
this  can  be  a time  saver.  In  practice  we  have  found  this  to 
work  quite  well  for  a range  of  velocities,  but  not  for  the 
time  increments,  which  depend  on  the  window  increments  for 
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resolution.  Averaging  of  the  data  reduces  the  random  com- 
ponents, but  must  be  done  with  discretion.  It  is  very  seldom 
that  the  return  signals  do  not  vary  to  some  extent  from  shot 
to  shot,  even  in  the  best  of  conditions.  Nature  never  quite 
follows  our  assumtion  of  flat,  laterally  homogeneous  layers 
of  sediments,  and  we  rapidly  begin  to  lose  information  if 
we  average  very  many  data  sets.  Again  experience  with  real 
data  provides  the  final  answer,  and  we  have  had  some  of  our 
best  results  without  averaging  over  shots,  and  summing  over 
frequency  only  in  the  final  stage  of  the  estimator.  After 
examining  the  statistics  of  the  estimators  in  the  next  chapter, 
we  investigate  the  results  of  applying  the  estimators  to 
real  data  in  Chapter  6. 
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Chapter  5 Statistics  of  the  Estimators. 

Introduction 

In  this  chapter  we  examine  the  statistical  distributions 
of  both  the  conventional  and  MLM  velocity/depth  spectral 
estimators.  We  calculate  the  bias  and  variance  of  both  forms 
of  the  estimator  for  Gaussian  input  data.  Simplified  results 
are  presented  for  the  special  case  of  independent  (between 
channels  and  between  observations)  noise.  With  the  aid  of 
a matrix  whitening  process  we  solve  the  estimator  forms  and 
their  statistics  for  the  case  of  a singular  (rank  1)  estimate 
of  the  covariance  matrix.  The  moments  for  the  MLM  estimate 
operating  on  a singular  covariance  matrix  are  shown  to  be 
a form  of  confluent  hypergeometric  function.  These  are  cal- 
culated and  compared  with  the  conventional  moments  using 
the  same  covariance  matrix  estimate.  The  MLM  is  shown  to 
improve  the  velocity/depth  spectral  estimate,  even  when 
employing  a singular  covariance  matrix. 

Throughout  this  chapter  we  characterize  the  data  (the 
output  from  the  FFT  operation)  as  a complex  valued  signal 
plus  complex  Gaussian  noise.  The  signal  is  considered  as 
unknown  but  constant,  and  the  noise  is  multi-variate  normal 
with  zero  mean  and  a covariance  matrix  £.  The  data  from 


IH 


i 


observation  "k"  is  denoted  by 


Xu  - 5 + M, 


where  S is  the  vector  of  signals  and  N is  a noise  vector. 
This  is  a common  assumption  in  geophysical  data  and  permits 
us  to  calculate  and  compare  the  statistics  of  the  two  forms 
of  the  estimators. 


Conventional  Estimator  Statistics 

We  begin  by  considering  the  statistics  of  the  conventional 


estimate.  We  use  matrix  notation  to  simplify  our  calculations. 


The  data  may  be  considered  an  NxL  matrix  formed  from  L obser- 
vations of  vectors  composed  of  the  N data  channels. 


T X - Tl 

X X.  X 
X X 


= U(f) 


The  estimated  covariance  matrix  is  formed  by 


R(f)  = U (fluff) 


5.3a 
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th  /\ 

where  the  ij  component  of  R is 

A L * 

R,.  * L 5.3b 

J W 

Referring  to  section  4-B,  the  averaging  inherent  in  this  form 
of  the  estimated  matrix  is  the  quadratic  averaging  in  column 
3 of  Figure  4.9.  This  is  the  form  generally  used  for  the 
estimate  of  the  covariance  matrix  of  a process  (Anderson  1958, 
Goodman  1963). 

We  write  the  data  as  a signal  plus  Gaussian  noise, 

= s}  + nifc  5.4 

where  is  distributed  as  N{£3,£),  = [s.jj  is  the  unknown 

but  constant  signal,  and  £ is  the  actual  covariance  matrix 
of  the  noise  process.  It  has  been  shown  (Goodman  1963, 

1965,  and  others)  that  the  estimated  covariance  matrix  j*  has 
a non-central  complex  Wishart  distribution.  This  is  a multi- 
variate generalization  of  the  non-central  complex  chi-square 
distribution. 

It  can  be  shown  that  for  a Wishart  distributed  matrix 
R and  a column  vector  of  constants  E,  the  quantity 
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Pc  - ¥ [ £+Be] 

is  distributed  as  a first  order  Wishart  with  L degrees  of 

x 

freedom,  which  is  equivalent  to  'X(L)  (Rao,  1965).  Specif- 
ically, the  distribution  is 

or*  'XV, \)  5.5a 


where  rrX 

E 

- 

5.5b 

and  Ais  a non 

-centrality  parameter 

given  by 

ibE^r 

5.5c 

If  we  look  at 

the  on-axis  response 

( Ej  = t ) an^  consider 

the  noise  to  be  uniform  and  independent  (Z.=G^I  ),  then 

we  have,  for  a 

simplified  case. 

__  a 

— Or 

N 

5.6a 

A 

__ 

5.6b 

The  characteristics  and  moments  of  the  non-central  chi-square 


97 


are  well  known.  The  mean  is 


e[R]  = 


and  the  variance  is 


<rS  = 0^(2.*?  A) 


For  our  simplified  case  we  have 


e[?]  = ov‘(j!r  + 1S) 


5.7a 


5.7b 


5.8a 


0-4/Ak  + 

UV  l,  N*  NO-*  J 


5.8b 


MLM  Estimator  Statistics 

We  begin  our  investigation  of  the  statistics  of  the 
MLM  with  results  derived  by  Capon  and  Goodman  (1970) . Using 
a relation  given  by  Rao  (1965) , we  can  derive  Capon  and 
Goodman's  result  in  a simple  fashion.  Rao  (1965)  gives 
the  following  result  for  a matrix  R distributed  as  W(L,£)  . 


The  quantity 


£+  £~'e 

E + R1  E 
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where  E is  any  fixed  vector,  is  distributed  as  7((L-N+1). 

The  numerator  is  a variance  term  which  remains  constant. 

We  rewrite  this  result  to  give 

I 

EL*  “ e'r-’e 

is  distributed  as 

i 1 <X*(L-N+l)  ; 5.9 

, i E E J 

This  result,  as  well  as  the  existence  of  the  Wishart  distri- 
bution density  function,  depends  on  L^N.  The  expressions  for 
the  mean  and  variance  are 

E[P«J  = [eVe]  (L-N+l+X)  5.10a 

= [e+£  e]  +4A)  5.10b 

where  'X  is  the  non-centrality  parameter  given  in  Eq.  5.5c. 

We  again  simplify  for  the  case  of  51*0^ J.  and  Sj  * S . The 
results  are 

r n / LrN*l  I 

E[P«u.]  = ” N 


L XtX  \ 

CTYX  ) 5.11a 
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F 

I 


E 


and 


01 


5.11b 


Comparing  Bq.s  5.8  and  5.11,  we  see  that  both  the  expected 
value  due  to  the  noise  power  and  the  variance  are  reduced 
by  the  MLM.  This  verifies  the  concept  that  the  MLM  is  a 
higher  resolution  estimator  and  does  not  respond  as  greatly 
to  incoherent  signals. 

We  note  in  applying  these  results  that  the  requirement 
that  L be  greater  than  N is  a problem.  In  most  of  our  appli- 
cations we  have  used  only  one  or  several  shots  or  observations 
in  forming  the  covariance  matrix.  In  order  to  examine  the 
statistics  of  these  cases,  we  propose  another  approach  which 
is  presented  in  the  next  section. 


MLM  for  a Singular  Covariance  Matrix 

The  general  results  of  Capon  and  Goodman  (1970)  are 
not  valid  for  the  case  where  L <N,  which  is  a region  we  are 
most  interested  in.  The  rank  of  the  estimated  covariance 
matrix  R is  L,  and  R does  not  exist  when  L <N.  We  have 


found  that  by  adding  a small  real  quantity  to  the  diagonal 
elements  of  6,  we  eliminate  the  singularity  of  the  inversion 
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and  can  invert  the  modified  matrix  even  when  its  original 
rank  was  unity.  This  operation  is  commonly  done  in  spectral 
analysis  techniques  and  is  described  as  whitening  the  matrix. 
The  effect  of  whitening  is  much  the  same  as  the  quadratic 
averaging  of  many  observations,  each  of  which  contain  some 
white  noise.  The  diagonal  terms  (which  are  all  zero  phase) 
are  enhanced  relative  to  the  off-diagonal  terms  (which  have 
non-zero  phases  that  vary  with  the  noise  components) . We 
note  that  the  diagonal  terms  are  spectral  components  and  the 
off-diagonal  terms  are  cross-spectral  components.  White 
noise  contributes  to  the  level  of  a spectral  estimate  with- 
out affecting  the  cross-spectral  level.  We  also  note  here 
that  the  linear  averaging  in  columns  1 and  2 of  Figure  4.9 
reduce  the  noise  by  increasing  the  number  of  observations, 
but  do  not  contribute  to  increasing  the  rank  of  the  matrix. 
The  estimated  matrix  following  the  linear  averaging  of  the 
terms  is  of  rank  1.  The  comments  in  the  next  section  on 
the  single  observation  case  are  also  applicable  to  this  case 
if  we  consider  a reduced  input  variance. 

Some  insight  into  the  singular  covariance  matrix  may 
be  gained  from  factoring  the  matrix  into  its  eigenvectors 
and  eigenvalues. 


r 


R 


WA  W 
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5.12 


TV  is  a diagonal  matrix  whose  diagonal  terms  are  "X,,  and  Wis 
a matrix  of  column  eigenvectors.  The  eigenvectors  are  ortho- 
normal so  that  wV=I  . R is  hermitian,  which  implies  that 

A 

the  Aj  are  real  and  non-negative.  The  inverse  of  R is  given 
by 


R 


w a"w+ 


5.13 


diagonal  matrix  whose  terms  are  (ty-  men  one  or 
of  the  Aj  are  zero,  R is  singular  and  the  inverse  is 
ill-defined. 

If  we  now  consider  the  modified  covariance  matrix 


R + pi 


5.14a 


we  get  a modified  eigenvalue  matrix. 

R = W A W + 01  5.14b 

or  R = W f A + 01)  W+  5.14c 


The  eigenvalues  of  the  modified  matrix  are 


102 


A;  = A,  + P 


5.15 


and  the  eigenvalues  of  the  inverse  are 


Ai+P 


5.16 


These  are  always  finite  for  0>O  and  the  inverse  is  no  longer 
ill-defined. 

The  stability  of  the  inversion  operation  is  determined 
by  a quantity  known  as  the  condition  number,  ^ (Householder 
1964) . The  condition  number  of  the  covariance  matrix  is  the 
ratio  of  its  largest  to  its  smallest  eigenvalue. 


^ = mr 


min  X 

j J 


5.17 


For  approaching  one,  the  inverse  operation  is  a well  posed 
problem  and  the  solution  involves  very  stable  calculations. 
For  % increasing,  the  computation  of  the  inverse  becomes 
more  and  more  unstable,  and  the  matrix  approaches  a singular 
condition.  By  whitening  the  matrix,  we  are  limiting  the 
value  of  The  largest  eigenvalue  is  bounded  by 


5.18 


103 


The  minimum  eigenvalue  in  the  modified  matrix  is  greater 

than  or  equal  to  £$.  Letting  ^ be  a function  of  Tr(R),  we 

control  the  stability  of  the  inversion  operation.  The 

optimum  maximum  value  of  ^is  determined  by  the  numerical 

stability  and  accuracy  of  the  computational  device  used 

for  the  calculations.  For  computers  with  a 7 significant 

4 

figure  accuracy,  a value  of  about  10  is  suggested.  The 
tradeoff  we  make  for  stabilizing  the  matrix  inversion  is 
one  of  distorting  the  matrix,  and  ultimately  distorting  the 
final  estimate.  Our  results  which  follow  indicate  that  this 
distortion  is  toward  the  conventional  form  of  the  estimator; 
hence  we  have  a valid,  if  slightly  less  than  optimum,  estimate. 
This  bias  toward  the  conventional  estimate  is  intuitively 
correct  in  that  as  we  increase  the  level  of  white  noise  in 
the  signal  field,  the  weighting  coeffecients  approach  uniform 
and  the  optimum  beam  pattern  approaches  the  conventional. 

We  conclude  that,  with  the  whitening,  it  is  possible  to 
employ  the  MLM  estimator  on  a singular  covariance  matrix. 

The  estimate  suffers  from  not  having  reduced  the  random 
components,  but  the  adaptive  procedure  should  still  produce 
a higher  resolution  estimate  than  the  conventional. 


4 
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Single  Observation  MLM  Estimator  - Exact  Solution 

A simple  case  that  we  can  solve  exactly  and  generally 
is  for  the  inversion  of  the  whitened  covariance  matrix  from 
one  observation  of  the  data  set.  The  estimated  covariance 
matrix  becomes 


T 

B = IT  + PI 


5.19 


where  ^ is  the  power  of  the  added  white  noise.  Gererally 

. _4  _2 

we  let  p be  in  the  range  of  10  to  10  times  the-  trace 
of  TTt  From  Graybill  (1969)  we  have 


V- 


YY 


5.20 


Solving  for  both  estimators  in  terms  of  the  vector  components. 


we  obtain 


P = — + -4  T t E*Y-  Y*E- 
K N N3-  L f-  E‘  T*  j J 


5.21 


X ' * 

ft  

NB  + N Zpi*  - it 

1*1  i */  j •(  J J 


5.22 
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We  make  the  following  substitutions 


V = 


£A, 


5.23a 


Ha, -a 


5.23b 


where 


A,  - e,  X 


5.23c 


A = jf  T.  A, 


5.23d 


We  note  that  V is  the  square  of  the  sample  mean  of  the 
steered  data  and  0 is  the  sample  variance.  Substituting 
these  into  Eqs . 5.21  and  5.22,  we  obtain 


JJL_ 


5.24 


5.25 


These  results  are  derived  in  Appendix  II.  We  use  these  ex- 
pressions for  the  estimators  in  the  remainder  of  this  dis- 


cussion. 


For  |3  increasing,  we  note  that  the  MLM  estimate  asymp- 
totically approaches  the  conventional  estimate.  The  two 
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estimates  also  converge  for  the  case  when  the  variance  of  the 
steered  samples  approaches  zero.  This  corresponds  to  the 
correctly  steered  estimate  of  a signal  without  noise. 

Statistics  of  Single  Observation  Estimators. 

In  order  to  determine  the  statistics  of  the  single 
observation  case,  we  again  consider  multi-variate  non-zero 
mean  Gaussian  data  as  the  input  to  the  estimators. 

T ~ N(£,£) 


5.26a 


A 


N ( [SiE*]  [EftSEj]  ) 


5.26b 


Constraining  the  noise  to  be  identically  distributed  and 
independent,  and  the  phase  corrected  signals  to  be  identical 
for  each  channel,  we  have 

A — 5.27 


Deriving  the  distributions  of  (f/ and  0 is  straight  forward, 
and  we  take  our  results  from  Papoulis  (1965) . is  a first 
order  non-central  chi-square  process.  The  probability  density 
function  is 


I 
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(±i£l 

HTy/N 


tjv_  i/V 

e^+ 

i 


U(f) 


5.28 


0 is  central  chi-square  with  density  function 


L- 


ttJ 


_ Ng> 
x<r; 


^(Wnw 


e 1 e T U(B) 


5.29 


We  also  have  that  ^ and  0 are  independent.  Changing  variables 


$ - 1 * F® 

© - f(f-0 

we  obtain 


and 


*LIA 


= ± + HL 
N f 


5.33 


Since  (fJ  and  ^ are  independent,  we  can  solve  for  the  distri- 


bution of  the  quotient  % 


<V  1 


in  a straight  forward,  but 


algebraically  complicated  manner  (see  Appendix  III) 


This  expression  is  not  reducible  in  terms  of  closed  form 
functions.  Because  and  J are  independent,  we  can  solve 
for  the  mean  and  variance  of  the  estimator  most  simply  by 


returning  to  the  expressions  for  the  distribution  functions 
of  ^ and  j . Calculating  the  first  and  second  moments  of 
and  -L  , we  have 

E[^J  = + f* 

L ‘ rv  5.3 

eM  = g + f 4 5.3 

E[f)=  (-kf  ViW.k) 


Efj^  " |?)  5.35d 

^(a ,b,x)  is  the  second  form  of  Kummer's  function,  a type 
of  confluent  hypergeometric  function.  (Abramowitz  and  Stegun, 
1965).  It  is  convenient  when  looking  at  the  mean  to  drop 
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1 

the  constant  % and  consider  only  the  mean  of  ^ for  the 
MLM  estimator  and  the  mean  of  for  the  conventional.  The 
mean  of  the  MLM  estimator  is  given  by 

E[ftJ  ‘ E[V]E[f]  5.36 

The  variance  of  the  estimator  is  given  by 

Va*[£J  = E[fi  ~ 5.37 

Using  the  relation 
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At(a  , b,x)  is  the  first  form  of  Kummer's  function,  which  may 
be  determined  form  tabulated  values  of  calculated  from  the 
series  definition  (see  Abramowitz  and  Stegun,  1965). 

The  expected  values  of  the  MLM  estimator  versus  the 
noise  variance  are  plotted  in  Figure  5.1.  As  the  noise  in 
the  data  increases  - as  the  signal  structure  deviates  from 
the  form  decreed  by  the  steering  vectors,  the  output  of 
the  estimator  drops  off  sharply.  This  is  an  indicator  of  the 
increased  resolution  that  we  find  with  the  adaptive  estimator, 
and  how  the  resolution  is  controlled  by  the  level  of  $ . As 
f}  is  decreased,  the  estimator  permits  smaller  and  smaller 
deviations  from  the  desired  signal  structure  in  order  to 
maintain  the  same  level  of  output  (i.e.,  the  resolution  is 
increased) . 

The  variance  of  the  single  observation  MLM  estimator 
is  plotted  in  Figure  5.2  along  with  the  variance  of  the 
corresponding  conventional  estimator.  Since  the  adaptive 
estimator  is  dependent  on  the  signal  field,  we  have  plotted 
the  variance  versus  the  variance  of  the  input  data  for  several 
values  of  signal  strength  and  for  different  values  of  additive 
white  noise.  For  small  values  of  0^.  the  variance  of 
is  highly  dependent  on  the  strength  of  the  signal.  As  the 
noise  variance  increases,  the  dependency  on  the  level  of 


VARIANCE  OF  ESTIMATE 
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additive  white  noise  becomes  the  dominant  factor.  For  com- 
parison purposes  we  also  have  plotted  the  variance  of  the 
conventional  estimator  for  the  same  signal  levels. 

Summary 

In  our  probabilistic  models  we  have  over  simplified 
the  actual  processes,  and  the  accuracy  of  the  models  suffers 
accordingly.  Without  these  simplifications,  however,  the 
problem  becomes  intractable.  There  are  too  many  dimensions 
to  allow  reasonable  interpretations  to  be  made.  Our  most 
vulnerable  simplifications  are  the  restriction  on  the  noise 
to  be  independent  and  the  restriction  to  looking  only  at 
correctly  steered  signals.  The  noise  in  seismic  reflection 
data  is  generally  highly  colored  and  propagating.  As  we 
scan  the  steering  vectors  versus  Tq  and  C,  what  was  the 
desired  signal  in  one  instance  is  a strong  part  of  the  noise 
field  in  the  next.  But  the  simplifications  do  allow  us  to 
identify  some  of  the  major  characteristics  of  the  estimator 
and  make  some  simple  comparisons.  The  results  from  the 
singular  matrix  case  are  applicable  in  a more  general  sense 
if  we  consider  unwanted  signals  as  non-random  components 
which  reduce  the  degrees  of  freedom  of  the  matrix.  Off-axis 


signals  increase  the  sample  variance  of  the  steered  data 
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and  increase  the  value  of  0 in  Eq.  5.27.  The  resulting 
estimator  shows  a much  higher  resolution  in  depth  and  velocity 
than  does  the  conventional  estimator,  which  is  one  of  the 
strong  points  of  the  adaptive  processor. 

As  we  saw  in  the  beginning  of  this  chapter,  the  multi- 
observation matrix  provides  its  own  whitening  while  reducing 
the  random  components  within  the  data.  Considering  the 
constraints  used  to  develop  the  estimator,  this  whitening 
should  be  optimum;  i.e.,  the  white  noise  level,  which  we 
have  seen  to  be  an  indicator  of  the  resolution,  is  a direct 
function  of  the  noise  variance  of  the  data.  Noisier  data 
calls  for  a broader  resolution.  The  single  observation 
MLM  estimator  is  sub-optimum  in  two  senses.  First,  there 
is  no  way  for  it  to  distinguish  between  signal  and  noise  - 
there  is  no  averaging  to  reduce  the  random  components. 

Second,  the  white  noise  level  is  externally  adjusted  and 
is  not  directly  related  to  the  noise  variance  of  the  data. 

In  spite  of  these  short  comings,  the  MLM  single  observation 
estimate  is  a large  improvement  over  the  conventional  single 
observation  estimate.  The  single  observation  MLM  estimator 
still  has  a greatly  increased  resolution  and  a reduced  var- 
iance. We  simply  need  to  be  judicious  in  our  estimate  of 
the  noise  level  of  the  data  and  our  choice  of  the  parameter  . 
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If  we  consider  noise  as  errors  in  the  data,  then  there  is  a 
tradeoff  between  the  mean  square  value  of  the  errors  which 
we  assume  to  be  present  in  the  data  and  the  resolution  that 
we  call  for.  We  cannot  afford  a high  resolution  if  there 
are  significant  random  or  systematic  errors  present. 

For  the  MLM  estimators  in  general,  we  find  that  we  have 
an  increased  resolution  and  a reduced  variance.  For  the 
single  observation,  and  any  case  where  the  random  components 
have  not  been  completely  averaged  out,  the  resolution  and 
variance  are  gained  at  the  expense  of  an  increased  bias. 

This  bias,  however,  which  is  a decreased  signal  level,  is 
a direct  contribution  to  the  increased  resolution  of  the 
estimate.  If  used  cautiously,  it  is  to  be  desired  rather 
than  avoided. 
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Chapter  6 Experimental  Results  and  Conclusions 

Introduction 

The  experimental  development  of  the  MLM  estimator  in  our 
study  has  been  through  a step  by  step  procedure.  We  began  by 
doing  studies  of  the  estimator  response  to  an  ideal  covariance 
matrix  as  in  Figure  2.9.  Following  this  investigation  we 
created  synthetic  tapes  in  which  we  could  control  the  charac- 
teristics of  the  reflectors  and  in  which  we  could  be  certain 
that  the  data  matched  the  travel  time  model.  Studies  with 
this  data  brought  out  the  requirements  for  careful  windowing 
and  demonstrated  the  increased  velocity  resolution.  Studies 
with  real  data  have  substantiated  the  results  and  conclusions 
of  the  synthetic  data  studies,  and  have  also  identified  some 
interesting  points  that  were  not  accounted  for  in  the  synthetic 
data.  These  include  the  ability  to  resolve  overlapping  sig- 
nals and  multiples,  and  the  strong  frequency  dependence  of 
certain  returns. 

Figures  6.1  and  6.2  give  typical  results  from  the  ideal 
covariance  matrix  studies.  Figure  6.1  is  the  conventional 
estimate  and  6.2  is  the  MLM  estimate  of  a four  reflector  case. 
The  reflectors  are  indicated  on  the  plots  by  a — j— . These 


spectra  were  calculated  entirely  from  a single  covariance 
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matrix  containing  all  the  reflectors.  The  ambiguous  strip  and 
the  higher  resolution  at  shallow  depths  that  we  predicted  from 
the  ambiguity  functions  are  evident.  The  effect  of  holographic 
focusing  is  indicated  by  the  resolution  in  both  depth  and  velo- 
city of  the  shallowest  reflector.  In  comparing  the  relative 
resolution  of  the  two  estimators,  the  MIM  processor  appears  to 
extend  the  focusing  range  of  the  array,  as  well  as  providing 
a higher  resolution  in  the  normal  range  of  the  conventional 
processor. 

The  next  step  in  the  development  of  our  study  was  to 
estimate  the  covariance  matrix  from  ideal  data.  A common 
ground  point  gather  from  a 12  channel  synthetic  tape  is  given 
in  Figure  6.3.  There  are  8 reflectors  whose  parameters  are 
given  in  Table  6.1.  This  data  does  not  contain  any  noise 
wavefronts  other  than  those  introduced  by  a filtered  random 
noise  generator  applied  to  each  channel.  Nor  does  it  contain 
multiples  or  refracted  arrivals.  All  the  delay  times  follow 
the  RMS  travel  time  model.  The  reflection  signatures  are 
damped  sinusoids.  The  conventional  and  MLM  spectra  of  this 
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Table  6.1 

Reflectors  in  Synthetic 

Data  Tape. 

(s«cl 

C (m/e) 

Primory 

Frequency  (Hz) 

.20 

1490. 

34. 

o 

00 

1840. 

30. 

1.30 

2260. 

27. 

2.10 

3050. 

24. 

2.50 

3230. 

22. 

2.70 

3490. 

20. 

3.45 

4120. 

18. 

3.60 

4430. 

15. 
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data  are  given  in  Figures  6.4  and  6.5.  These  estimates  are 
summed  over  frequency  and  include  the  frequency  band  from 
19  to  35  Hz.  An  interesting  point  of  comparison  is  the  flat- 
ness of  the  background  in  the  adaptive  estimate;  the  absence 
of  much  of  the  fine  structure  that  is  present  in  the  conven- 
tional estimate.  The  noise  is  white  Gaussian  and  we  are 
observing  the  reduced  variance  predicted  in  Chapter  5.  The 
resolution  along  the  time  axis  is  comparable  for  both  estimates, 
as  we  would  expect  from  Figures  2.8  and  2.9.  There  is  a 
notable  increase  in  the  velocity  resolution  in  the  MLM  esti- 
mate, particularly  for  the  deeper  reflectors.  The  spectra 
for  6 channels  of  this  same  data  ( the  even  numbered  channels, 
giving  the  same  array  length)  are  given  in  Figures  6.6  and  6.7. 
The  results  of  the  adaptive  procedure  applied  to  the  6 channel 
array  are  not  quite  as  good  as  in  the  12  channel  estimate, 
but  continue  to  be  greatly  improved  over  the  conventional 
estimate. 

Following  the  synthetic  data  studies,  we  turned  our 

1 

We  note  that  most  of  our  spectra  plots  are  contoured  at  6 dB 
intervals.  These  were  plotted  before  we  normalized  the  gain 
factors  between  the  conventional  and  adaptive  estimation  pro- 
grams, so  we  can  only  say  that  the  levels  are  arbitrary.  All 
the  contours  above  an  arbitrary  level  are  shaded  to  aid  in 
interpretation.  The  shading  contour  was  chosen  in  each  case 
to  highlight  the  peaks. 
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Figure  6 . 5 


Data  Adaptive  Velocity/Depth  Spectrum  of 
12  Channel  Synthetic  Data.  19-35  Hz. 

6 dB  contour  levels. 
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attention  to  actual  field  data.  The  data  is  from  the  WHOI 
multi-channel  system  employing  six  channels  at  150  meter 
spacings.  The  shot  points  are  spaced  at  37.5  meter  intervals, 
so  we  would  expect  a fairly  high  correlation  between  adjacent 
gathers.  The  next  four  figures  give  the  conventional  and 

I 3 

adaptive  spectra  of  two  consecutive  common  ground  point 
gathers.  Figures  6.8  and  6.9  give  the  conventional  estimates, 
and  Figures  6.10  and  6.11  give  the  adaptive  estimates.  We 
see  the  same  reduction  in  sidelobe  energy  and  flatness  of 
spectrum  that  we  observed  in  the  synthetic  data.  Note  the 
reflector  at  0.85  seconds  and  1800.  meters  per  second  that  is 
virtually  lost  in  the  energy  from  the  shallow  refracted  and 
direct  arrivals  in  the  conventional  spectra.  The  higher 
resolution  of  the  MLM  allows  it  to  discriminate  between  the 
direct  and  refracted  returns  and  those  returns  fitting  the 
RMS  travel  time  model.  Comparing  the  reflector  at  1 . 7 seconds 
and  2400  meters  per  second  in  Figures  6.9  and  6.11,  the 
adaptive  estimate  distinguishes  between  the  reflector  and  a 
slower  multiple,  while  the  conventional  estimate  smears  them 
together.  Finally,  the  general  stationarity  of  the  data 
between  adjacent  gathers  is  an  indication  that  we  can  further 


improve  the  estimate  by  averaging  the  covariance  matrices. 


4000 


Figure  6.8  Conventional  Velocity/Depth  Spectrum  of 
6 Channel  Georges  Bank  Data.  19-35  Hz. 
Shotpoint  300.  6 dB  contour  levels. 
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Figure  6.9  Conventional  Velocity /Depth  Spectrum  of 
6 Channel  Georges  Bank  Data.  19-35  Hz. 
Shotpoint  301.  6 dB  contour  levels. 
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Figure  6.10  Data  Adaptive  Velocity/Depth  Spectrum  of 
6 Channel  Georges  Bank  Data.  19-35  Hz. 
Shotpoint  300.  6 dB  contour  levels. 
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Another  example  of  the  improvements  gained  by  use  of  the 
adaptive  estimator  is  given  in  Figures  6.12  and  6.13.  This  is 
another  WHOI  6 channel  data  set  from  Georges  Bank.  We  again 
observe  the  same  relative  improvements  of  the  MLM  over  the 
conventional  estimate.  Note  that  Figure  6.13  is  plotted  at 
3 dB  increments  to  bring  out  the  structure  of  the  spectrum. 

The  energy  from  the  direct  and  shallow  refracted  arrivals  in 
the  conventional  estimate  is  greatly  attenuated  in  the  MLM 
spectrum.  Two  reflectors  at  1.45  and  2.25  seconds  and  2700 
meters  per  second  are  greatly  enhanced  relative  to  slower 
velocity  returns  in  the  adaptive  estimate.  The  velocity 
resolution  of  the  reflector  at  3.35  seconds  and  3100  meters 
per  second  is  significantly  better  in  the  adaptive  estimate. 

This  same  data  also  gives  us  a good  example  of  the  infor- 
mation partitioning  as  a function  of  frequency.  The  spectra 
in  the  previous  six  figures  (6.8  - 6.13)  have  all  been  aver- 
aged over  a 19  to  35  Hz  frequency  band.  If  we  look  at  the 
estimates  for  each  frequency  component  for  shotpoint  1020 
(Figure  6.13),  we  find  that  the  different  travel  paths  are 
highly  frequency  selective.  Monochromatic  MIM  estimates  for 
this  shotpoint  are  given  in  Figures  6.14  through  6.18.  Most 
of  the  multiple  energy  is  in  the  higher  frequencies  of  the 
band,  and  most  of  the  penetrating  primary  energy  is  in  the 
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Figure  6.13  Data  Adaptive  velocity/depth  spectrum  of 
6 channel  Georges  Bank  data.  19-35  Hz. 
Shotpoint  1020.  3 dB  contour  levels. 
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Figure  6.14  Adaptive  Velocity/Depth  Spectrum  of  6 Channel 
Data.  Frequency  Breakdown  of  Shotpoint  1020. 
19  Hz.  6 dB  contour  levels. 


Figure  6.16  Adaptive  Velocity/Depth  Spectrum  of  6 Channel 
Data.  Frequency  Breakdown  of  Shotpoint  1020. 
27.  Hz . 6 dB  contour  levels. 
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lower  frequencies.  This  is  indicative  of  the  filtering  done 
by  the  travel  path  medium,  and  points  out  the  usefulness  of 
working  in  the  frequency  domain  for  both  the  conventional  and 
adaptive  spectra  estimates. 

We  can  also  use  this  data  to  demonstrate  the  effects  of 
averaging  over  shots.  The  estimated  spectra  (MIM)  for  shot- 
points  1021  and  1022  are  given  in  Figures  6.19  and  6.20.  A 
spectral  estimate  using  shotpoints  1020  - 1022  in  forming  the 
covariance  matrix  is  given  in  Figure  6.21.  Some  of  the 
reflectors  (1.25  seconds)  are  improved,  while  some  (2.3  and 
3.35  seconds) are  degraded  over  the  better  of  the  single 
observation  estimates. 

| 

Conclusions 

We  have  seen  how  the  Maximum  Likelihood  Method,  when 
applied  to  velocity/depth  spectra  estimation,  gives  improved 
resolution  of  reflector  parameters.  The  resolution  in  depth 
is  determined  by  the  windowing  of  the  data,  which  is  identical 
for  both  the  conventional  and  MLM  processors.  The  resolution 
in  velocity  is  determined  primarily  by  the  coherent  power 
estimate,  and  here  the  MLM  introduces  significant  improvements. 
The  improvement  is  dependent  on  the  additive  noise  field,  but 


is  substantial  for  the  synthetic  and  field  data  that  we 
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Figure  6.20  Data  Adaptive  Velocity/Depth  Spectrum  of 
6 Channel  Georges  Bank  Data.  19-35  Hz. 
Shotpoint  1022.  6 dB  contour  levels. 
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examined.  The  increased  velocity  resolution  is  of  particular 
significance  for  deep  returns  where  the  array  is  approaching 
the  limit  of  its  near  field  geometry.  Here  the  MLM  can  extend 
the  practical  operating  range  of  the  array. 

We  have  identified  two  aspects  of  applying  the  estimator 
which  are  critical  to  its  successful  implementation.  The 
windowing  is  an  important,  although  subtle,  aspect  of  the 
estimation  procedure.  Its  dual  role  of  transform  window 
and  arrival  time  detector  creates  a tradeoff  between  time  and 
frequency  resolution.  The  expected  form  of  the  signature 
and  its  frequency  spectrum  must  be  considered  in  resolving 
this  tradeoff.  In  addition  to  the  windowing,  another  critical 
area  is  the  averaging  of  data.  Here  again  we  find  tradeoffs, 
this  time  between  the  stability  of  the  estimate  and  the 
ability  to  resolve  reflectors  that  may  be  present  in  only 
one  or  few  shots  and/or  frequencies. 

In  examining  the  statistics  of  the  estimates,  we  have 
extended  previous  results  to  include  our  application,  as  well 
as  derived  the  statistics  for  a new  case;  the  single  obser- 
vation MLM  estimate.  A by-product  of  this  is  a new  expression 
for  the  MLM  estimator  which  promises  to  be  a much  faster 
algorithm  to  implement  (Eq.  5.25).  The  MLM  has  a reduced 
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variance  which  is  obtained  at  the  expense  of  an  increased 
bias.  This  bias  can  be  controlled  somewhat  through  the  matrix 
whitening  parameter.  It  is  a direct  consequence  of  the  higher 
resolution  of  the  adaptive  estimator,  and  is  not  a severe 
problem  if  the  whitening  parameter  is  adjusted  to  suit  the 
data . 

In  concluding,  the  advantages  of  applying  the  MLM  esti- 
mator are  primarily  in  the  resolution  of  velocity,  both  when 
distinguishing  between  multiples  and  primaries  and  at  the 
limit  of  the  array's  operating  range.  The  method  allows 
the  use  of  shorter  and  sparser  arrays  to  obtain  results 
equivalent  to  a conventional  analysis,  and  hence  reduces 
operating  and  processing  costs.  On  the  other  hand,  it  resolves 
reflectors  that  would  not  be  resolvable  with  the  conventional 
analysis  for  a given  data  set. 


A 
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Appendix  I Development  of  a simplified  expression  for  the 
mean  square  error  of  the  Fourier  transform  due 
to  windowing  of  the  time  series. 

The  error  is  given  by  Equation  4.5 

® jafrPt 

m = + „(Dw(t-r,)]  e dt  4.! 


If  we  square  this  and  take  the  expected  value,  we  have 


CD  CO 

1 j-fft.-U 

+ e [n(tj  n(tj J w(t(-rjw(trrw)  J e 


AI.l 


noting  that  the  cross  terms  drop  out  because  we  have  assumed 
^(t)  and  n(t)  independent.  We  can  rewrite  the  expected  value 
terms  as  autocorrelation  functions. 


E|>(t = o(t,-TA)  R.fc-tJaftt-Z)  ai.2, 


E [n(ii)n(tJJ  = Rn(i.-tJ 


AI.2b 


We  can  write  the  autocorrelation  function  as  the  inverse 
transform  of  the  frequency  spectrum. 
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R ft.- 1.) 


e 


Al.3 


Substituting  this  into  Equation  AI.l,  we  have 

|E(«r  = 5xWpt,U,a(t,-'Ci)aftrirj  » 


+ 


, , -j 

dpSnw  U Jk"(t,-Z)"(tx-z)e 


j 


j 


AI.4 


The  integrations  in  t,  and  t factor  and  are  complex  conju- 

1 2 

gates  of  each  other. 

/ / , 

|e(|)|*  = U'S.M  jdto(t-T,){w(i-Tw)-lJe 


r 

( -j11 r(f-v>)£ 

dt  w(t,-rw)e 

AI.  5 


If  we  now  assume  that  Sx  and  S are  approximately  constant 


n 


for  (f-V)  small;  i.e.,  within  the  bandwidths  of  the  squared 
quantities,  we  can  take  them  out  of  the  integral  and  we  get 
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(()  f = S*  (f)  dt  °(t-D  ( w(tx)-ij  e 


Letting 


+ 5„(f)  jdv  |dt  w(t-r„)  e' 
g(t)  = a(t-%j{w(t-r,j-l) 


-j  ITT  (f-v)t 


AI.6 


Al.7a 


we  have 


h (t)  - w(t-rw) 


G(f-v)  = ( j(t)  e 


AI.7b 


AI.8a 


/ -jiTr(f-vJt 

H (f-W  = j h (t)  6.  c/6  AI.8b 


Substituting  this  into  Equation  AI.6,  we  then  have 


E(f) 


AI.9 


From  Parseval’s  theorem  we  have 


G(f-v) 


G(») 


= ( dy  NM 


AI .10a 


Al.lOb 


Substituting  AI.10  into  AI.9,  we  finally  obtain 


5 JO  (dt(°(t-Z)[w(t-rJ-l}) 

+ sn(f)  (ji(wii-rj)1 

* 
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Appendix  II  Derivation  of  simplified  expressions  for  esti- 
mators using  one  observation  of  the  data  set. 

The  two  estimators  are 


Pc  = 


All . 1 


All  .2 


with 


B = YT  + PI 


All. 3 


e *E  = N 


All. 4 


We  have 


IT 


All. 5 


Multiplying  out  the  conventional  estimator  and  gathering  the 
terms  in  summations,  we  have 


Pc  = > Z E. 


. N N $ * 


« w * 


= £ * ill  fE, 


All. 6 
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The  MLM  estimator  becomes 


This  gives  us 


Q i N N J 

P = + All  AA 

K‘  N N Si  £ 1 J 


All.  9 


P = 

MLM 


f>*  + P 5 a,  a; 

Np  + N 21  Mi 


- ££  A A* 

*■'  J-*  J 


All. 10 


We  can  further  simplify  PMTM  with  the  following  two  identities. 

MlxM 


All .11a 


E A A = I Mi 


E Z A, A*  = |EAt 


All. lib 


We  also  make  use  of  a simple  theorem. 


Theorem  I 


Z x,1  = E(Xi-x)1  + (Zk) 

wticre  X = jy  21  Xj 


All. 12a 


All. 12b 


Proof : 


.2.  N , z 


r(x,-s)  - r(x,  -zxx,  + 


x1) 


=.  f x,‘  - 2*1*1  + NX 
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Appendix  ill  Derivation  of  the  probability  density  function 
of  z = (V/f  ) . 

Both  y and  J are  constrained  to  be  positive.  We  are  there- 
fore limited  to  the  first  quadrant  in  the  plane.  We 

first  calculate  the  cummulative  distribution  function  F 


AIII  .1 


Differentiating  with  respect  to  2a,  we  obtain  the  density 
function 


AIII. 2 
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The 

and 


joint  density  function 


is  simply  the  product  of 


c 


rrx 


x 


r 1L& 

~ <rJ- 


<rY‘ 


UiHUM 


aiii. 3 

Substituting  this  into  Eq.  Am. 2,  we  then  obtain  an  integral 
expression  for  the  density  function. 


£(*)  = n ** 

j 


CO 


z 


AIII. 4 


This  does  not  appear  to  be  solvable  in  closed  form.  For 
the  case  of  ^«0  it  reduces  to  a form  of  confluent  hyper- 


geometric function. 
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Glossary  of  Notations  and  Symbols 


Notations 


N 


Complex  conjugate. 

Complex  conjugate  transpose. 

Transpose. 

is  distributed  as. 

Vector  or  matrix  quantity. 
Estimated  quantity. 


Vector  or  Matrix  with  elements  xi 


t 
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Symbols 


alt) 

A 

A 

P 

f 

R) 

C,  C; 

C 

d 

D 

AT 

E, 

E(  f) 
E L-] 

fe 

ff 


Envelope  for  stochastic  signal  model. 

Array  element  gains  for  adaptive  processor  (2). 
Steered  data  vector  (5). 

Coefficients  for  series  expansion  of  travel 
time  (1 ) . 

Matrix  whitening  parameter  - quantity  added  to 
diagonal  elements 

Uniform  signal  magnitude  for  simplified  statistics 
Complete  gamma  function. 

til 

Seismic  velocity  within  (i  ) layer. 

RMS  velocity. 

Array  spacing  (3). 

Depth  of  first  layer  (1). 

Difference  in  position  of  window  and  position 
of  center  of  desired  signal  (4) . 

Steering  phase  for  i*"*1  channel. 

Steering  vector  of  phase  shifts  used  to  focus 
the  array. 

Error  in  transform  due  to  windowing  and  to  noise  (4 
Expected  value  of  (.) 

RMS  bias  error  of  transform  due  to  windowing. 
Probability  density  function  for  0 . 

Probability  density  function  for  ^ . 


a 
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i 

0 

r 

k 

Kw 

X 

L 

A 

Ai 

A 

A\ 

n(t) 

n'j 

N 

Mk 

N(S,£) 

&n*0,WPc 


Reducing  coefficient  in  adaptive  estimator  (5). 
Sample  variance  of  steered  data  (5). 

Identity  matrix. 

Wavenumber  vector. 

Window  factor. 

Matrix  condition  number. 

Number  of  observations  used  in  forming  the 
covariance  matrix. 

Ray  parameter  (1).  Non-centrality  parameter  (5). 

Eigenvalues  of  R. 

Eigenvalue  matrix  of  j*. 

Half  width  of  data  window. 

First  form  of  Kummer ' s function. 

Noise  in  data. 

th  th 

Noise  from  i channel  and  j observation. 

Number  of  channels. 

th 

Noise  vector  from  k observation. 

Multi-variate  complex  Gaussian-normal  distri- 
bution. 

Conventional  velocity/depth  estimator. 

Maximum  Likelihood  Method  velocity /depth 
estimator. 


<P(Kek0-) 

<P(Z,cJt;c,) 
<LkcxlT c.) 
Ulti,  X\Q,  tC 
%(kA) 

V 

Rff) 

Aftef),  R(f) 

fd.r.clfj.c) 

Mi) 

J(  f) 


5(f) 

5x(f) 

5n(f) 

SMh) 


Angle  with  vertical  of  wave  vector  in  i 
layer  (1). 

Monochromatic  plane  wave  ambiguity  function  (3). 
Monochromatic  velocity/time  ambiguity  function. 
Wideband  velocity/time  ambiguity  function. 
Chi-squared  distribution  function. 

Plane  wave  steering  function. 

Square  of  the  sample  mean  of  steered  data  (5). 

Covariance  matrix  of  data  field. 

Estimated  local  covariance  matrix  for  data 
windowed  by  T,  C. 

Normalized  response  of  velocity  depth  array 
(conventional)  to  frequency  bandwidth. 

Reflected  signal  in  data. 

Fourier  transform  of  signal -4(t). 

Signal  from  channel  j. 

Vector  of  signals  in  the  data. 

Frequency  spectrum  of  seismic  source  signal  (3) . 
Frequency  spectrum  of  process  x(t). 

Frequency  spectrum  of  process  n(t). 

MLM  wavenumber  estimator. 

Covariance  of  noise  process  n(t),  channels  i and  j 
Distribution  scaling  factor. 
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Variance  of  conventional  estimator. 

Variance  of  MLM  estimator. 

Uniform  noise  variance  of  data  for  simplified 
statistics . 

Covariance  matrix  of  noise  process  N. 

Travel  time  (one-way)  through  i layer  (1). 

til 

Travel  time  thickness  of  i layer  (1) . 

Two-way  acoustic  travel  time  from  source  to 
(j  n)  receiver. 

Normal  incidence  two-way  travel  time. 

Trace  of  ( . ) . 

Half  width  of  signal  envelope  a(t). 

Delay  of  signal 

Delay  of  data  window  w(t). 

Unit  step  function. 

Data  matrix. 

Second  form  of  Xummer ' s function. 

Data  window  function. 

Wishart  distribution  function. 

Eigenvector  matrix  of  R. 

Wideband  stationary  process  of  stochastic 
signal  model. 

Source  to  ( j ) receiver  distance. 

Data  from  channel  i. 


Frequency  domain  representation  of  signal 
from  channel  i (observation  j). 


Y(f),  Yj 


Vector  of  channels  of  frequency  domain  repre- 
sentations (observation  j). 


Ytr.c-.f) 


Frequency  domain  representation  of  data  in 
windows  positioned  according  to  T,  C. 


Z 


Non-constant  part  of  adaptive  estimator. 


I 


Unitary  vector  (all  ones) . 
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